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THE GRAVITY STIFFNESS OF A SUSPENSION- 
BRIDGE CABLE 


By A. G. PUGSLEY (University Engineering Laboratories, Bristol) 


[Received 2 October 1951] 


SUMMARY 

This paper discusses the stiffness of a heavy inextensible cable in terms of the 
work done against gravity when the cable is loaded, and examines the relations 
between this energy treatment and the conventional ‘deflexion theory’ in common 
use. The latter is often presented in a form that appears to imply that the gravity 
stiffness of a cable is negligible; this is seen to be misleading and to result from 
neglect of a term in the expression for zero extension and of the discontinuities that 
arise with concentrated loads. General expressions for the gravity stiffness are 
given appropriate to applied loads small compared with the cable weight, and the 
gradual reduction in the accuracy of these expressions as the applied loads are 
increased is examined. 


1. THe mechanics of a suspension-bridge structure involves the inter- 
action of the two major components, the suspension cables and the bridge- 
deck girders. An early theory, developed by Rankine (1), was based upon 
the assumptions that the dead load—the weight of the whole bridge 
structure—was carried by the cables and that the live load—any vehicles, 
ete., passing across the bridge—was insufficient to cause any deflexion 
of the cables from their natural configuration appropriate to the dead load. 
The bridge-deck girders were thought of as being stiff enough to spread 
the live load effectively over the whole span of the cables. Rankine’s 
theory is still used for the approximate discussion of suspension-bridge 
problems, but the design of such bridges, particularly in the U.S.A., is now 
largely based on the so-called ‘deflexion theory’ of the suspension bridge, 
due originally to Melan (2). In this theory some allowance is made for the 
small deflexions that arise from the passage of the live load. The essentially 
non-linear response of the cable to such loads is realized, and some correction 
for it attempted. 

From the mathematical standpoint this deflexion theory, in its con- 
ventional form, has been very conveniently stated by von Karman and 
Biot (3), and their notation, apart from the adoption of v for the vertical 
deflexion of the cable in place of their w, will be employed. This change is 
made to bring the equations nearer those discussed in recent English 
papers, as by Atkinson and Southwell (4). The conventional theory notes 
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that a cable with a horizontal component of tension H carrying a. tot 


al 
vertical loading of q¢ per unit length (horizontal) will adopt a form de. 


fined by 
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” dat — —% () 


where x and y are horizontal and vertical coordinates measured from one 
end of the cable as in Fig. 1, and that the same cable, if called upon to carry 








nein ~ 








Pre. 1. 


an amount p, (per unit length) of the live loading p, will develop a vertical} 
deflexion v given by the equation 
d? 
(H +-h)— (y+v) -—Pr; (2) 
dax2 
where / is the extra horizontal component of the cable tension due to p,. 
The remainder of the live load, that is (p—p,), must be carried by the 
beam (assumed to be simply supported at its ends), and if this is uniform, 
the corresponding equation of equilibrium will be 
d4y 
El )— pj. (3) 
dx! } } 1 
Combining (1), (2), and (3), we have as the fundamental equation of this 
deflexion theory: 
d4y dv h 


dat (Hh) Ts p~4 H’ “ 


El 

In this equation there are two unknowns, the function v and the parameter 
h. The further condition necessary for the solution is that the change in 
length of the cable, due to the deflexion v, must be consistent with the 
increase fh in its tension and its extensibility. In many cases the cable 
is relatively inextensible and a first approximation corresponding to a 
+ Horizontal deflexions, which are neglected in the conventional theory, have been dis- 
cussed by Atkinson and Southwell (4). 
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GRAVITY STIFFNESS OF A SUSPENSION-BRIDGE CABLE 387 
condition of zero extension is used. This condition is expressed approxi- 
mately by the relation r 
| wg dx Q), (5) 
0 
where L is the span of the cable. 

2. From the foregoing statement of the current theory in its simplest 
form it will be seen that the problem has been resolved into the simultaneous 
solution of the differential equation (4) and the integral equation (5), and 
that emphasis is thus placed on the problem of estimating the increase of 
tension h. An alternative viewpoint has been discussed elsewhere (5), (6) 
by the author. In this the stiffnesst of the cable is seen largely to arise from 
its own weight and the dead load upon it, and the interaction between 
cables and deck girders is regarded as the interplay of the gravity stiffness 
essentially non-linear) of a cable and the elastic stiffness (linear) of the 
beam. This emphasis on a gravity stiffness brings to the fore a type of 
stiffness that is not generally familiar, but it helps to present the suspension- 
bridge problem in clear physical terms and in a way that brings out the 
nature of the non-linearity present. 

It is natural first to look at the problem of such a cable in terms of the 
effect of a single concentrated load, as illustrated in Fig. 1. If the dead 
weight gL of the cable (including its suspended dead loads) is large com- 
pared with the applied load (P), then the deflexions (vertical and hori- 
yontal) of any point on the cable will vary linearly with P; but as P becomes 
comparable with gL these deflexions will become large enough to alter the 
general geometry of the cable, and the deflexions v, and wu, will tend towards 
limiting values (if the cable is inextensible) set by the straight line (dotted) 
configuration indicated. Thus the ‘gravity stiffness’ of the cable, measured 
by ¢P ev,, will increase with P as this becomes of the same order as qL. 
The work done by P as it is applied and moves through v, (= }Pv, for the 
linear case) is stored. not as elastic strain energy in the cable, but as gravity 
potential energy; in other words, application of the ‘live’ load P raises the 
centre of gravity of the dead load system associated with the cable. Thus, 
for the linear case we may write 

L. 
1Pr, | vq dx. (6) 
0 

But we have seen in (5) that the conventional approach treats the integral 

in (6) as approximating to zero, as though the gravity stiffness of the second 


approach were negligible and the cable in a state of neutral equilibrium! 


Here considered as a stiffness against vertical loads. 
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This paper is hereafter devoted partly to the resolution of this difficulty. 
and partly to the estimation of the non-linear response of the cable to 
vertical loads. 

3. It is clearly desirable to examine first more critically the use of the 
expression (5) for an inextensible cable. In the presentation given by von 
Karman and Biot it is noted that the initial total length is given by 


=f b+(ayi : 


and thence, by replacing y by (y+-v) and expanding in a Taylor series, the 
variation AS is found to be 


ms (8) 


j (dy/dx)(dv/da) d. 
{1+ (dy dx)?! 


neglecting higher terms in v. Thence, by integration by parts, noting that 
v = 0 at the two limits of the integral, and neglecting the departure of the 
denominator of (8) from unity, it is found that 
L L 
: ad?) o£ 
AS a [ Y dx : vg dx, (9) 
J dx H., 
0 0 
by substitution from (1). Equation (5) is just the particular form of 
(9) when AS is zero. 
From a review of the approximations made in this argument it is appar- 
ent that the paradox noted has been unaffected by the replacement of 
-(dy/dx)*\} in (8) by unity; the explanation must lie at least in part in 








: eae : , a | 
the neglect of higher terms in the Taylor series for S. If we include a 


further term we have 


L L 
: - (dy\ (dv 1 ¢ (dv\? 
S Be dx 4 lx. 10 
” | Falta "Ts | (7) Pe 
0 


The second term here is, of course, the same as the change in length, due 


to v, of a straight member from 2 = 0 tox = L. 
Again, the replacement of the first of these two integrals in (10) by 
L 
(1/H) | vq dx as in (9) depends upon a process of integration by parts which 
0 


is strictly legitimate only when dy/dx and dv/dx are continuous functions | 


of x between the limits x = 0 and x = L. In general, this will not be so 
unless the live loading p, on the cable is also such a continuous function, 4 
point which seems to have been ignored by most investigators. It 1 
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obviously not of much importance when the live loading is distributed 
over a considerable proportion of the span and when the presence of a 
stiffening girder ensures that the loading p, on the cable is still further 
distributed, but in the case of an isolated cable loaded by a single con- 
centrated force, the adoption of an identity between 


I L 


. . 


[ (SY)(52) ae and 7 vg dx 


ey 
0 0 


is erroneous and may have serious effects. 
Thus, for our case of a loaded cable under a small single load P, we have 


as the two governing equations the energy equation (6), i.e. 
L 
4 Pv, | vq da, (11) 
0 
ind, for an inextensible cable, the further equation 


L 


FNC) ans f(a =o a9 


0 0 


A special case arises when P is applied at the centre of the cable span. 


We can then concentrate on the half span and note that 
I L bL +L 

1 dv d . a * dy 1 f¢ 
(- J | dx Jy “7 y dx - y vdx = — | vqdx. 

J da da d r 0 3 ( x . C x H J 


0 0 0 





(13) 
The integration by parts is legitimate in this case because of continuity 
between the limits, at each of which either dy/dz = 0 or v = 0. For this 
special case only, therefore, can the expression (12) be written 


I bL 


ef (fe 0 


0 0 


thence. from (11 


we have the direct relation 
+L 1 9 
1p av “ad. 7 m 
1Pv, = H | (7) (15) 


There is thus a real grav ity stiffness associated with a loaded cable that, 
ior values of P small compared with qZ, is defined by (11) and (12) in the 
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general case and by (15) for the case of central loading. If we are concerned. 
as in the analysis of section 1, with a distributed live loading p, but on an 
isolated cable, then (11) will be replaced by 


3 - 
vp dx | vg dx, (16) 
0 0 
L 
it being assumed that the live load | p dx is small compared with the dead 
0 
L 
load | q dx. 
0 

Possible reasons why equations such as (11), (12), and (15), and the limi- 
tations of equations (9) and (10), have not hitherto been noted, are that 
the conventional deflexion theory envisaged distributed rather than con- 
centrated loads, and that modern investigators have only returned to the 
consideration of such loads for a complete cable-plus-girder system, and 
not for a cable alone. It is interesting to note that Timoshenko’s bridge 
analysis (7), upon which much modern work has been based, treated v as 
expressible across the whole span by the first few terms of a Fourier series, 
and so automatically envisaged continuity. 

4. So far as work in this country is concerned, it is necessary to go back 
to the middle of the nineteenth century to find any relevant systematic 
published work on the deflexion properties of a loaded cable. The problem 
of a uniformly loaded cable, hanging initially in parabolic form, and then 
loaded by a single concentrated load P, is easy to solve formally, but 
intractable algebraically when an explicit expression for its deflexion due 
to P is required. Whilst seeking to produce such an expression of an 
approximate nature, the writer came upon the work of an unnamed author 
that was published in a little known engineering journal in 1862 (8). This 
paper examined, inter alia, the behaviour of a cable under the conditions 
discussed above. The dead loading q is taken to be uniform across the 
span L and a small concentrated load P is applied at a general point along 
the cable. Coordinates x and y are measured from an origin at the lowest 
point of the cable, of dip d, and the vertical load P is applied at x = rL. 
The cable is treated as inextensible and the only approximation made is 
that this condition can be expressed by (10) with the second term omitted. 
We will return to this point later. 

It is shown that if the equation of the initial parabola is 
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the parameter ™, becomes, for the two new parabolic ares into which the 
application of P causes the cable to deflect, m,+-m,, where 


2. eT (18) 
sf 2qL 


(Consistent with this, the increase of horizontal tension becomes 


or | " 
h Als 7 )i—4. (19) 
= 
rT 6 
| } 


4 4 Lt 








Of more immediate interest, the deflexions of the point where the load P 


s applied are 
Pd (1+-12r?)(1—4r?) (20) 
2qL 1+(3P/2qL)(1—4r?)’ 7 


‘ SP a" r(1— 167%) = (21) 
qL? | (SP 2g L)(1 47°) 
For a central load, 7 0, and these expressions become 
ha . (22) 


0 SgL 143(PiqL)’ 


and when, in addition, P is very small compared with qL, 


Pd 


23 
2q a (59) 


It is of interest to study this last result in relation to equation (15) of the 
special case of our energy discussion. The 1862 paper referred to gives, 
lor that case, 

(l+n)? LZ? (x+4nL) , 2 


d _d . —_— a, ; (24) 

m,(1+8n) 4 m,(1+3n)  m, 
wheren = P/qL. This represents a parabolic variation of v with x as shown 
n Fig. 2. The point of zero vertical deflexion occurs at x = 3L, and the 
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aL +L 
areas | vdx and | v dx to either side of this point are numerically equal 
0 tL 
+L 
but of opposite sign. The net area | v dx is thus zero, in accordance with 
0 
(13), and the assumed equation for inextensibility. 


Each of these component areas has the value 


dL{ n m 
27 \1+ 3n)’ (29) 


and we may note that in (14), for this case of central loading with q constant, 


the first integral iL 


q[ » da, 
H 


0 


in which H is qL*/8d, is the difference of two equal terms, each amounting to 


8d? n 98 
27L\1+3n)' °°) 


This is to be compared with the value of the second integral in (14), which 
for the same case is, from (24), 


1L 
] | dv ¢ d? | ail i (21) 
2 J \dx L\1+3n 

0 


It will be seen that for small values of n (= P/qL) the terms (26) are much 
larger than (27); but that as m increases the two approach each other in 
value. + 

We thus see that for a heavy cable under a relatively small live load the 


assumption of L 
| dy t ——" 
J \dx}\dn 


0 
as the condition of inextensibility involves the difference of two nearly 
equal quantities, each much larger than the second term in (12) that is 
commonly omitted. This near-balance of two large quantities accounts 
for the tendency of the cable to behave somewhat like a see-saw under 
quarter-point loads, though the cable is, of course, always stable and the 
gravity stiffness always positive. Treatments ignoring the second term 
in (12) are likely to give a good approximation only when the live loading is 
small compared with the dead loading, and is well distributed along the 
span. Moreover, from the energy treatment for a concentrated central 


+ Expressions (26) and (27) become equal when n = 0-53, but neither is accurate to so 
large a value of n. 
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load, the inclusion of the second term in (12) for such a case is evidently 
essential. In that case, for » very small, (27) gives 


I 


| (ase = 8) 


0 
Thence, from (15), which applies exactly under these conditions, 


Pd 
2qL 


(] 
0 


4s in (23). 











0-1 fl UM = ind - --(23) 
+ ai 
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5. We have in the equations (20), (21), and (22) a first step towards 
expressing the non-linearity of the gravity stiffness of a loaded cable, 
and it is of interest therefore to consider (22) in relation to the obvious 
limit that must arise as P/qZ becomes large. In this case, the cable will 
take up a triangular form, and if it is inextensible it is easy to show that vp 
must become asymptotic to 4d. The form (22), on the other hand, tends to 

las P/gL + oc, so that as the deflexion becomes large equation(22), whilst 
satisfactorily indicating an increasing stiffness as P/qL increases, does not 
indicate a large enough increase. In general, expressions (20), (21), and 

22), from the nature of the assumptions made, would not be expected to 
be accurate for values of P/qL greater than about 0-1 or 0-2, and this is 
confirmed by experiment. The 1862 paper pursues this matter to obtain 
wccurate, but indirect, expressions for the cable deflexions when P/qL 
exceeds these limits. Fig. 3 is based upon its results, and shows, for a 
centrally loaded inextensible cable with a dip equal to one-tenth of its 
span, the accurate relation between P and v as compared with the linear 
telation (23) and the approximate non-linear relation (22). It will be seen 
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that the approximate relation follows the accurate one closely up to 
values of n = P/qL of at least 0-1. The general equations (19), (20), (21), 
and (24) thus provide a useful approximate representation of the non-linear 
response of a cable to vertical loading. 
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THE FLOW OF WATER UNDER A SLUICE-GATE 


By A. M. BINNIE (Trinity College, Cambridge) 
[Received 7 August 1951] 


SUMMARY 


The analysis of the two-dimensional free motion of a perfect liquid under a sluice- 
ein an open horizontal channel is simplified by introducing the Froude number, F, 


the stream. The relation is obtained between F, and F,, the values of F on the 
pstream and downstream sides of the sluice-gate. If, as is usual, F, < 1, then 


1. Dimensionless expressions for the depths of the stream, the discharge, and 
horizontal force on the gate are found in terms of F,. Experimental and theoreti- 
nvestigations into the necessary gate opening are correlated by plotting them 
» base of F,, and Rayleigh’s result for flow under pressure through an aperture 


s shown to hold good up to a remarkably large value of F,. An explanation is put 


forward for the absence of waves on both sides of the sluice-gate, from which it is 


luced that a Scott Russell wave of depression is impossible. 


1. Introduction 

MANY investigators have studied the two-dimensional motion of water 
inder a sluice-gate in an open horizontal channel. Consequently plentiful 
nformation has appeared concerning the discharge when the upstream 
onditions and the sluice opening are specified. In their published form, 
however, these results apply only to the precise circumstances in which 
they were obtained, and no correlation exists between them. On the 
theoretical side, Southwell and Vaisey (1), in the course of their compre- 
hensive attack on free-surface problems, showed how the velocity distribu- 
tionin the channel can be determined by relaxation methods. As anexample 
they gave the complete results for one particular case. 

The purpose of this paper is to show how greatly the problem may be 
simplified by the introduction of non-dimensional variables—a procedure 
which was used by Rouse (2) and by Black and Mediratta (3) in their 
treatment of the hydraulic jump. All the necessary calculations have been 
performed once and for all with the exception of the determination of the 
gate opening, and here the experimental and theoretical results mentioned 
i the preceding paragraph have been harmonized to provide the answer 


to this question also. 
2. Possible types of flow 

We consider two-dimensional motion in a straight, horizontal, and rect- 
ngular channel of constant width unity, which conveys a uniform stream 
01 frictionless liquid, the discharge and the total head (measured above the 
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bottom of the channel) being denoted by Q and H. Rouse (4 and 2) and 

Southwell and Vaisey (1) showed that the depth, d, of the stream can haye 

two values. This conclusion is easily proved, for the equation of continuity 

1S Q vd. (2.1) 

v being the uniform horizontal velocity, and Bernoulli’s equation is 

ye 

H—d4 _, 
29 





125 














oS) 


pe 
Fic. 1. Roots of the equation (2+ F?)? = KF?, 


where g is the acceleration due to gravity. The elimination of v leads to 


~, 
®Ha+” — 0, (2.3) 
— 2g 
which is a cubic equation ind having in general two positive roots lessthanH. 
[t is, however. more informative to introduce the dimensionless Froude 


number F, defined by i 
7 


oe aes (2.4) 
(gd)* 
In place of (2.3) we then obtain another cubic equation 
(2+ F?)8 = KF?, (2.5) 
in which the non-dimensional quantity K is given by 
, 89H? 
Kk = “9 — (2.6) 
* 


The nature of the roots of the cubic equation is seen in Fig. 1, where the 
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curve z = (2+ F?)® is plotted on a base of F*, together with one of the 
straight lines 2 KF? cutting the curve at A and B. Since negative values 


of F2 are inadmissible, that part of the diagram is omitted. The limiting 


position of the straight line z — KF? is the tangent which touches the curve 
it the point F l,2 27. The roots of (2.5) are given by the intersections 


4and B; therefore there are in general two possible values of F’, one greater 
ind the other less than unity. The velocity V of a long gravity wave of 
small amplitude moving on the surface of liquid of depth D is 


} (gD); (2.7) 


nd on comparing (2.4) with (2.7) we see that at F 1 the stream possesses 

this velocity. There is an urgent need for suitable words to describe the 
conditions when F is greater and less than unity. The adjectives ‘shooting’ 
ind ‘streaming’ have occasionally been employed, but they are far from 
iid. Certain classical scholars, accustomed to coin words for scientists, 
havesuggested ‘superundal’ and ‘subundal’, and these terms will be adopted 
throughout this paper; they are analogous to ‘supersonic’ and ‘subsonic’, 
which are used in accounts of the motion of gases. 


Fig. 1 shows that the analysis yields no result if K < 27, i.e. if 
Q > (39H). (2.8) 


The failure can be explained by supposing the flow to originate in an infinite 
reservoir of depth H, the channel being regarded as an extreme form of 
broad-crested weir. This motion is a particular example of critical flow. 
The general solution was shown by Binnie (5) to be a stream of such a depth 
that the velocity is equal to that of a long gravity wave on the surface; in 
these circumstances the discharge is the greatest that can be obtained with 
the given head H. In this instance the depth is 2H/3 and the velocity is 
29H /3)' in conformity with (2.2) and (2.7). Hence the maximum discharge 
is (89H*/27)*, and a greater value of Q may not be specified. 


3. The « fect of the sluice-gate 

When a sluice-gate is lowered into a superundal stream and fixed in 
position, a violent commotion accompanied by splashes occurs in front of 
the gate. (The proviso ‘fixed in position’ is important, because Southwell 
and Vaisey (1) showed that a Scott Russell solitary wave can be formed on 
the upstream side of a sluice-gate so held that its bottom edge is above the 
level of the undisturbed stream.) Attention will therefore be confined to the 
usual case of an approaching stream possessing the lower value of F, say F,, 
which is yielded by (2.5). To facilitate the calculation of F, from specified 
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values of Q and H, Fig. 2 has been prepared in which K is displayed on, 


base of F,. Logarithmic plotting has been adopted in order to give pronj- 


nence to low values of F, which are those commonly encountered jy 


practice, but natural values are also shown. It should be noticed that, wher 


F? is very small in comparison with 2, (2.5) reduces to 


R\1 
F, = ft 
= (5) 
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Fic. 2. Variation of K, F,, and d,/d, with F,. 


We shall be concerned only with the state in which the discharge is 
unobstructed, and a side elevation of the channel and sluice-gate under 
these conditions is shown in Fig. 3, in which symbols with suffix | apply to 
the upstream side and those with suffix 2 to the downstream side. The 
surface of the stream rises appreciably in the neighbourhood of the sluice- | 
gate, and a stagnation-point C, where the velocity is zero, occurs at a height 
H above the channel bottom. Downstream the surface falls away until, 
theoretically at infinity, it again becomes horizontal, and a uniform distti- 
bution of velocity is restored. Here the stream has the larger value of F, 
say F,, which is given by (2.5), and necessarily it moves with a superundal 


velocity. Thus the sluice-gate may be regarded as a device which trans- 


forms a subundal into a superundal stream. 
The variation of F, with F,, obtained from 


(2+F3? (247) 


2 L 


F F 


is also displayed in Fig. 2. Again, an approximation similar to (3.1) may be 


made, for when F3 is large compared with 2 


F, = Ki. 
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d on Since (2.1) and (2.4) may now be written 
pron V V} d, Uede, (3.4) 
red ys ry 
F a" R= —-, (3.5) 
. Whe) ; (gd,) r (gd.)' 
VU; 
———fee H 
| g T 
v2 
= Y Y 
aM Fic. 3. Side elevation of stream near sluice-gate. 
, 2 
J e 


ve find, on eliminating v, and vy, that 


. d, i 2 2 @ 
C . ad 3.4 
? d, (7) ‘ (9.6) 


[he values of the ratio d,/d, have been added to Fig. 2, and with the aid 





fthis diagram the properties of the issuing stream may readily be deter- 
ined when Q and H are given. 
| 4. The formation of waves 
‘8° 8 | The results shown on the previous pages throw light on the question 
ee: whether trains of standing waves can occur on either side of the sluice- 
.* | gate. In a slightly different form this matter was raised by Rayleigh (6), 
: ‘ whose results were reviewed by Lamb (7, sections 242-5). It was examined 
prs iain by Southwell and Vaisey (1), who expected to discover waves as their 
| culations proceeded; in fact, none was found on either side. This outcome 
“ ‘ {their work requires scrutiny, and since possible waves may have infinitely 
7a small or finite amplitude, both types must be considered. 
ia The well-known analysis of the former class takes the profile of the waves 
or | tobesimple harmonic. On water of uniform depth D their velocity depends 
o ipon the wave-length A and, in accordance with (2.7), has the maximum 
lue (gD)! when A is great. Hence these waves cannot keep up their 
position on the stream below the sluice-gate where F > 1. 
(3.2) | To deal with the question of finite waves on the downstream side, use 
, will be made of Allen’s theorem which was expounded by Southwell and 
vy be Vaisey. This theorem states that the surface can nowhere fall below the 


evel d,. The stream now being recognized as superundal, the theorem can 


le expressed more generally as follows: No stationary wave possessing a 
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trough can exist on a superundal stream. In this form the theorem ca) 
readily be proved by a reductio ad absurdum argument. We suppose that 
as shown in Fig. 4, a finite trough exists on a superundal stream of depth 
d, and velocity v,. At the trough the depth and surface velocity are denote; 
by d and v’. Since the stream-lines in the trough are concave upwards and 
possess the same Bernoulli constant, the velocity over the cross-section L}) 








Fic. 4. Wave trough on superundal stream. 


is greatest at the free surface L (Fig. 4). Hence the equation of continuity 
leads to 


; -Vede . i‘ f 
v x2 in which X > 2. (4,] 
d 

mM y'2— ve ve r9 d > ve 79 9 

Thus ee ——] — (X2x?— 1) (4,2 
2g 2g\ & 2g 
: , : d, 
where we write, for brevity, x 7° (4.3 
€ 


But Bernoulli’s equation shows that 


vie—ve og. ( bot *) (44 
2q x x 


therefore the right-hand sides of (4.2) and (4.4) are equal. Then, since 
v5/(gd.) > 1, we find, after some reduction, that 


3a—2 > X72, (4.5 
The straight line z = 3x—2 an:| the cubic z = 2° are shown in Fig. 5, the 
irrelevant part where x < 1 being excluded. Since z = X2z° lies above 


the cubic, the condition (4.5) cannot be satisfied, and we conclude that 
the original assumption is false. Thus a train of finite waves possessing 
crests and troughs cannot occur on the downstream side of the sluice-gate. 

It may be remarked in passing that the extended theorem explains why 
Scott Russell’s solitary wave can be found only as a wave of elevation. The 
analysis given by Lamb (7, section 252) shows that the velocity of this wave 
is restricted to a small range of superundal values. The wave being regarded 


as stationary on a uniform stream, the lowest velocity of the stream is given 
by & = 1 when the wave is of zero height. At the other extreme, McCowai 
(8) proved that for the highest wave F = 1-25. Neither theory nor exper 
ment has discovered a similar wave of depression. 
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For the upstream side different arguments are required. Close to the 
sluice-gate the velocity over a vertical cross-section of the stream is least 
at the surface; and since Ff, < 1 this disturbance caused by the sluice-gate 
extends for an infinite distance upstream. Thus if a train of small simple- 
harmonic waves were superposed on the surface, the discharge through the 
elevated positions would be less than that through the depressions no longer 


8 
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ae 


Fic. 5. Proof of extended theorem. 


occupied by liquid; and in order to maintain the discharge, the mean level 
would rise against gravity. We deduce that the waveless surface is stable. 

A train of finite stationary waves is precluded by the fact that if it existed 
it would, because F;, 1, extend into the reservoir at infinity and disturb 
the conditions there. This matter is worth pursuing farther, and it will be 
demonstrated by a different method of reasoning that no known type of 
stationary wave train can occur. To a second approximation the theory of 
finite progressive waves on liquid of finite depth was obtained by Stokes (9) 
in Cartesian coordinates x and y. Many years afterwards (10) he extended 
the work to the third approximation but expressed it in the form of equations 
giving « and y in terms of the velocity potential and the stream function. 
This approximation in terms of x and y was placed on record by Bowden (11). 
Ifwe are content with the second approximation we find that the discharges 
across vertical planes under the crests and troughs are equal. But as Stokes 
pointed out in his later paper, the motion disclosed in the third approxima- 
tion is radically different, for there is now a drift, greatest at the surface, 
of the particles in the direction of wave propagation. Consequently, as was 
stressed again by Bagnold (12), the velocity of propagation is indeterminate 


and various definitions may be used. Thus the artifice of bringing the waves 
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to rest by superposing even a non-uniform constant velocity may not be 
employed, and it seems that this type of wave train cannot maintain itself 
undistorted in front of the sluice-gate. In their study of stationary finite 
waves Southwell and Vaisey postulated that the undisturbed streay 
possesses uniform velocity; thus these waves too are excluded. 

No permanent waves were found in the course of experiments in a channe| 
6 in. wide and 3 in. deep, supplied with a uniform stream by means of the 
orifice method described by Binnie (13). A train of progressive waves could 
readily be generated by external means in the length of channel between 
the sluice-gate and the orifice; but when the stimulus was removed, the 
train died away in the same manner as if it had been produced in a rectangu- 
lar tank of stationary water. It was noticed also that the theoretical rise 
in level to the stagnation point on the sluice-gate did not occur, this effect 
being very marked at large values of F,. Whether the sluice-gate was 
vertical or inclined forward or backward, in front of it a sharply defined 
zone of eddying water was formed, the surface of which moved in the 
reverse direction to the main stream. This exemplifies the well-known 
instability of viscous fluid motion against an adverse pressure gradient. 

The foregoing theoretical remarks throw some light on the discrepancy 
reported by Southwell and Vaisey. Working with an undisturbed stream 
having a constant value of F = 0-426, they obtained the profiles of station- 
ary waves of various heights. Thus they were able to plot wave-length 
against crest height; and by means of a small and reasonable extrapolation, 
they deduced the wave-length of an infinitely small wave. This considerably 
exceeded the wave-length given by the usual formula 


; ga 2nd 
y2 — tanh ; (4.6 
2a A 
which in non-dimensional form may be written 
. A 2rd re 
F2 — tanh : (4.7) 
27d A 


[t now seems clear that the finite waves which they investigated must have 
a shape differing somewhat from Stokes’s waves, but it is surprising that, 
when the amplitudes of both are reduced to zero, the wave-lengths should 
differ by as much as 6 per cent. 

Sufficient information is available for limits to be placed upon the size 
of stationary waves on a stream of velocity v, and depth d,. The crest 
height A, measured above the undisturbed depth d,, may not rise above 
the stagnation level; hence the maximum value of A is given by 
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With the aid of (3.5) this expression becomes 
Ow (4.9) 
d, 2 
A restriction on the trough depth B, again measured from the level of the 
undisturbed stream, is set by Allen’s theorem, from which we find that the 


sreatest value of B is given by 


B st (4.10) 
d, d, 
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Fic. 6. Maximum size of stationary wave. 


nd this can be evaluated by means of (3.6). The results are shown in Fig. 6 
ona base of F,, A d, being plotted upwards and B/d, downwards. Southwell 
0-426. At so 


small a value of F, Allen’s theorem is ineffective, presumably because the 


ind Vaisey’s result is indicated by the vertical line at F, 


wave-length is small, hence the stream-lines at the trough possess consider- 
able curvature and v’ in Fig. 4 is greatly in excess of the mean velocity over 
the cross-section of the trough. At larger values of F,, and hence of the wave- 
ength, it is possible that, although A becomes greater, the limit prescribed 


by the theorem is more nearly approached. 


5. The resultant horizontal force on the gate 

A force P (Fig. 3) is required to maintain the gate in horizontal equili- 
brium, and its magnitude may be determined by means of the principle 
oi momentum. The block of water, enclosed between vertical cross-sec- 


tions situated far away on each side of the sluice-gate, is exposed to three 
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horizontal forces, namely P and the hydrostatic forces on the two cross. 
sections. On equating the net force to the rate of change of momentum, we 
obtain Po] 
4(d2—d2) —_— = — (v3d,—v?d,), (5.1) 
p g 
where p denotes the density of the liquid. When this equation is combined 
with (3.4) and (3.5), the result is 








P = al 14 F2 » s+ F3 | 5.3 
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Fic. 7. Horizontal force on sluice-gate. 





Its correctness at both extremes can readily be checked. When the sluice- 
gate is lifted nearly clear of the stream and F, > F,, P > Oasit should. Ifthe 
sluice-gate is lowered to obstruct the flow almost completely, F, > 0 and 
F, > «©; then P + 4$pH?, which is the hydrostatic force on the shut gate. 
With values of F, obtained from (3.2), the dimensionless quantity P/(4pH?) 
calculated from (5.2) is plotted in Fig. 7 on a base of F,. 

This diagram is of little importance to a hydraulic engineer, but it is of 
interest because it can be used without modification to determine the drag 
of a planing surface having any two-dimensional form. The problem of 
planing surface advancing with velocity v, in water of depth d, is identical 
with that of a sluice-gate fixed in a stream approaching with velocity v, and 
depth d,. Up to the present no mention has been made of the shape of the 
sluice-gate, hence all the results obtained in the preceding pages hold good 
for the planing surface also, provided that the water is shallow. A well- 
known distinction is drawn between ‘long’ waves and surface waves in deep 
water, and for similar reasons a theory dependent upon Froude number 
annot be expected to apply without modification to a body only slightly 
immersed in deep water. 
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6, The opening of a vertical sluice-gate 

We now turn to the outstanding difficulty of determining the sluice 
opening S, and here the simple methods hitherto employed offer no assis- 
tance. Recourse must be had to experiment and, on the theoretical 
side, to relaxation methods. It is convenient to work with a ratio, such as 
§/H or the coefficient of contraction d,/S. The former is preferable because 
t is of greater practical importance, it is easier to measure, and more 
experimental values have been published. The stream approaching the 
sluice-gate is completely specified when v,, d,, and g are given, and a 
straightforward application of dimensional theory shows that S/H must 
a function of the Froude number only. Accordingly, in order to deter- 
mine this function, experimental and theoretical results will be plotted on 
: base of F,. 

Gibson’s experiments (14) were concerned with a sluice-gate 3 ft. wide 
nd with three openings, 0-1063, 0-2063, and 0-3083 ft., which he designated 
{, B,and C. From his curves of coefficients of discharge, the values of Q, 
H, and K were calculated for heads above the top of the orifice equal to 
1, 0:2, 0-5, and 1-0 ft. The corresponding values of F, were then found 
from (2.5). The twelve points thus obtained are shown in Fig. 8, which is 
livided into two parts, the results for S/H < 0-16 being plotted separately , 

Addison (15) used a model 6 cm. wide with three gate openings (a) 1-99, 
)) 2-98, (c) 6-00 em. and total heads varying from about 14 to 53 cm. The 
published diagram of his results is small; therefore a table of his actual 
observations was used, which he kindly provided in order that the highest 
possible accuracy might be obtained. For each of the three sets, four pairs 
f observations of Q and H were selected, uniformly spaced over the range 
ftotal head; and after analysis as explained in the previous paragraph they 
were added to Fig. 8. The theoretical result obtained by Southwell and 
Vaisey (1) has also been inserted. 

In addition, use may be made of the theoretical treatment due to Ray- 
eigh (16), who found by means of a conformal transformation that, for a 
stream emerging under pressure in the absence of gravity, 


d. 

2 — (2-611. (6.1) 
I . ° ° ° . . 

When the motion is due only to gravity, the corresponding case is given by 


SH +> 0,v, > 0, d, > H, F, + 0, d, > 0, v3/(2g) > H, F, > 0. We also see 
irom (3.1) and (3.3) that F, F3 > v8; thus (3.6) becomes 


(6.2) 
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Then, on combining (6.1) and (6.2) and identifying d, with H, we find that 
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F, are very small, is shown in Fig. 8. 


Rayleigh’s result is seen to hold good up to remarkably large values oi F,, 


Southwell and Vaisey’s point, which was obtained for F, 
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(0-426, is only 


just above the straight line, showing that even there the lack of uniformity 


of pressure over the upstream side of the sluice has a very small effect: it is 
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compensated by a slight increase in what may be termed the Rayleigh 
sjuice-opening. In general, the experimental results lie above the line by 
small amounts which increase slowly as F, becomes greater; and boundary- 
laver effects require the sluice to be raised a little distance additional to that 
necessary in frictionless flow. But for small values of F, Rayleigh’s result 
is very accurate, as was pointed out by Addison (15), who, however, did not 


proceed beyond (6. l .. 


REFERENCES 
1. R. V. SourHWELL and G. VaisEy, Phil. Trans. Roy. Soc. A, 240 (1946), 117. 
2. H. Rouse, Elementary Mechanics of Fluids (London: Chapman and Hall, 1946). 
3. J. Buack and O. P. Meprratta, Aero. Quart. 2 (1951), 227. 
4. H. Rouse, Fluid Mechanics for Hydraulic Engineers (McGraw-Hill, 1938). 
5. A. M. Brynte, Proc. Roy. Soc. A, 197 (1949), 545. 
6. Lorp RAYLEIGH, Proc. London Math. Soc. 15 (1883), 69 ; Collected Papers, 2, 258. 
7, H. Lams, Hydrodynamics (Cambridge, 1932). 
8. J. McCowan, Phil. Mag. (5) 38 (1894), 351. 
9, G. G. Sroxes, Trans. Camb. Phil. Soc. 8 (1847), 441; Collected Papers, 1 (1880), 
197. 
10. Collected Papers, 1 (1880), 314. 
11. K. F. Bowpen, Proc. Roy. Soc. A, 192 (1948), 408. 
12. R. A. BAGNOLD, J. Inst. Civ. Eng. 27 (1947), 447. 
13. A. M. Brynte, Aero. Research Committee, R. & M. No. 1887, 1940. 
14. A. H. Grsson, Proc. Inst. Civ. Eng. 207 (1920), 427. 
15. H. Appison, J. Inst. Civ. Eng. 8 (1938), 53. 
16. Lorp Ray.eicH, Phil. Mag. (5) 2 (1876), 441; Collected Papers, 1, 297. 














THE DECAY OF SHOCK WAVES 


By W. CHESTER (Dept. of Mathematics, The University, Bristol) 
[Received 25 September 1951] 


SUMMARY 

The interaction of the shock and rarefaction waves produced in a gas-filled tube, 
closed at one end by a moving piston, is discussed. The piston is moved impulsively 
to produce two waves of rarefaction separated by a shock. The decay of the latter, 
resulting from its interaction with the two plane waves, is described on the basis of 
an approximate treatment in which all terms of the third order in the shock strength 
are ignored. In particular, entropy variations are disregarded. It is shown that the 
strength of the shock decays inversely as the time. 

The application of the results to the more general problem, in which the piston 
experiences a series of impulsive changes in velocity in such a way that shock waves 
and centred rarefaction waves are formed alternately, is discussed. 

Finally the analogous problem involving the two-dimensional steady flow of a gas 
along an indented wall is considered. 


Part I. ONE-DIMENSIONAL UNSTEADY FLOW 


Introduction 

THE basic problem to be investigated is that shown diagrammatically in 
Fig. 1. A piston, originally at rest in a gas-filled tube, is moved impulsively 
to produce two plane waves separated by a shock wave. The subsequent 
behaviour of the shock wave is studied, assuming that the rarefaction waves 
are of sufficient strength to ensure that their interaction with the shock wave 
does not terminate within a finite time. 

A complete description of the subsequent gas-flow is a most difficult 
problem, not simply because of its non-linear character but, more particu- 
larly, because of the variations in entropy resulting from the propagation 
of a shock of non-uniform strength. However, if the shock strength is not 
too large, a second-order solution of the problem can be obtained neglecting 
the entropy variations, since these are of the third order in the shock strength 
(1, p. 156). Such a simplification has already been applied by Friedrichs (2) 
to the decay of a shock bordered on one side only by a rarefaction wave. 

The resulting theory should be a good approximation within a limited 
distance from the piston, but it is of interest to study the asymptotic be- 
haviour of the shock. The present theory, for example, shows that in the 
problem envisaged in Fig. 1 the strength of the shock decreases inversely 
as the distance from the piston. The validity of such a result may be 
questioned on the grounds that the neglected third-order variations in the 
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THE DECAY OF SHOCK WAVES 409 
flow quantities, produced by the interaction, accumulate to variations of 
a more significant magnitude when the time becomes indefinitely large. In 
an attempt to overcome such an objection, a critical examination of the 
validity of such asymptotic expansions has been made by Lighthill (3). His 
conclusions suggest that when the flow ahead of the shock is known exactly, 
«9 that conditions on the other side can be specified with an error of order 




















Fic. 1. 


the cube of the local shock strength, then the second-order theory is valid 
for all time in the sense that the error in the velocity of the shock is of the 
third order. Such is the case in the problem of Fig. 1, since any interactions 
which occur cannot affect the flow ahead of the shock. 

A further consequence of Friedrichs’s theory (2) is that, within the approxi- 
mation considered, the only effect of the interaction between the shock and 
rarefaction waves is a decrease in the strength of the shock, and the flow 
on either side of the latter continues to be described by the original simple 
waves until they are annihilated at the shock. Thus, for any impulsive 
motion of the piston in which shock and rarefaction waves are produced 
lternately, the subsequent behaviour of the shock waves may be deduced 
at once from the solution of our basic problem, since the simple wave regions 
remain unchanged by the interactions to the present order of approximation. 
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The example illustrated in Fig. 2 is a particular case, except that the firg 
shock will be propagated into a region of undisturbed gas. The behavioy 
of such a shock has already been discussed by Friedrichs (2). 

There is, however, an important qualification to the results describing the 
asymptotic behaviour of the shocks in this more general case. Condition 
ahead of each shock subsequent to the first are no longer known exactly; in 
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fact there is an error of order the cube of the maximum initial shock strength. 
In this case, Lighthill concludes that the present theory describes the final 
stages of decay correctly only toa first approximation. This, of course, does 
not imply that the analysis is no better than linear theory, which cannot 
deal with the problem at all since the strength of the rarefaction waves is 
not small, and some interesting conclusions can be obtained in spite of this 
restriction. 


Notation and basic formulae 
We consider in detail the following problem. A piston, originally at rest in 
a gas-filled tube, is given a succession of impulsive velocities —u,, —u , —Us 


where w, is positive and w, and uz are greater than u,. The initial retraction 
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hegins at time ¢ = t, when the piston is at a point x = 2, of the tube, and 
initiates a centred rarefaction wave which penetrates the undisturbed 
vas (Fig. 1). The origin of the space-time coordinates is taken to coincide 
with the point at which the piston velocity changes from —u, to —u,. Since 
i, > Uy, aShock wave will be sent out from the piston with uniform velocity, 
eventually overtaking and interacting with the rarefaction wave proceeding 
from (,, t,). Finally, a second rarefaction wave is initiated at the point (z,, t,) 
of the (x, t)-plane where the piston velocity changes from —u, to —u . The 
points at which the two rarefaction waves first meet the shock path are 
denoted by (x,,¢,) and (2,,¢,) respectively. For the present it is assumed 
that the wave proceeding from (,,t,) meets the shock path first, so that 
t <%g,t, <t,. A sufficient condition for this to hold is t, > —t,. 

It is convenient to quote at the outset some well-known results governing 
the behaviour of rarefaction waves and shock waves. 

Let a centred rarefaction wave radiate from the point (x,,t,) of the (a, t)- 
plane and let the particle and sonic velocities in the wave be denoted by u,, ¢, 


respectively. Then, inside the wave (see 2), 


r—2X, = w,(t—t,), (1) 
u, = (1—p2)(w,—¢9) (2) 
¢ Cotp?(w,—Co), (3) 
Ww, = U,+C, (4) 
p : 4 (5) 


where y is the adiabatic index, cy is the velocity of sound in the undisturbed 
fluid, and w, is a parameter which is constant along each ray of the wave. 
In particular, if the particle and sonic velocities in one of the uniform flow 
regions bounding the wave are —w,, c, respectively, then equations (2) and 
(3) may be replaced by 


U,+Uy (l—p?)(w,—w,), (6) 
C,—C, p(w, Ww), (7) 
where Wy C;— Uy. (8) 


Now let the suffixes r, s denote the conditions on the two sides of a shock 
wave having a velocity U, then (see 1, p. 149) 


and, correct to the second order, 
u,—u, = (1—p?)(w,—w,), (10) 
C.—C, = p?(w,—,), (11) 
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since, to this order of approximation, the transition relations through a 
shock agree with those for a simple wave of equal strength (1, p. 156). 

A combination of equations (9), (10), and (11) yields the following equa- 
tions, correct to the second order, 


9 
7 ; v.— U,—U,)* 
U = ou, pete 4 (tet) ” 


30 —p?) * 8e,(1— 








2 (WW) ++ = (w,—w,)?. (13) 
SC) 


We now return to our original problem. In the initial stages of the motion 
(¢ < t,) the velocity of the shock is constant, and the shock path is therefore 
represented in the (x, t)-plane by a straight line through the origin separating 
two regions of constant state in which the flow velocities are —u,, —u,, 
respectively. The point (,, ¢,) is thus given by the intersection of the shock 


path, namely st nate OE 


U;—Uy , (Uy—Us) 
rT 9) Be,(1— 2)?” 
with the tail of the rarefaction wave emanating from (x,,t,) which is, from 
equation (1), 


9 


where U Cc 





(14). 


Y—X, = (c,—Uu,)(t—t,), (15) 
since W, = Wy = C,— Uy, (16) 
along this path. 

With the help of the further relation, 


Lt, = —U,1,, 
, St - 
one can easily show that t= — r, (17) 
: a(4+-c) 
Uu Us ) 
where >. 8... (18) 
c,(1—p*) 


Interaction of the shock wave with the first rarefaction wave 
Between the points (a,,¢,), (a2, ¢,) the shock path will be modified by its 
encounter with the rarefaction wave emanating from (x,, ¢,). In general this 
interaction will affect the flow to the rear of the shock, but, as previously 
stated, the deviation from the uniform state already existing behind the 
shock will be a third-order quantity in the shock strength. Since the flow 
pattern ahead of the shock cannot be influenced by the interaction, the 
flow pattern on the two sides must remain unchanged to the second order. 
Furthermore, equation (1), which describes the flow ahead of the shock, 
may also be regarded as the equation of the shock path in which w, is a func- 
tion of ¢ and varies continuously along the path. Comparing this equation 
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with the expression for the gradient of the shock path deduced from equation 
(13), namely, dx, 42 ( )2 (19) 
4(w,-+-w,) + — (w.—w,)?, ¢ 
dt 7 " Sc,* ” . 
we obtain the following differential equation 
d ] , . 
mies t.)} 4(w,+we)+ = (w,—w,)*, (20) 
( < Cy 
in which w, is the independent variable and t,, w., c, are constants. The 
solution of (20) is easily found to be 


4ke 
Ww, Ws - Bs ey : (21) 
{(t—t,) /(t, —t,) }#—k 
where the constant of integration has been chosen to make w, = w, when 


t = t,, and 


Equation (1) now gives for the shock path 


Cc —t.)t(t—t.)3 
= #,+w,(t—t,) 4ke,(t, —t,)*(t—t,) ; 
; 1—k{(t,—t,)/(t—t,)}* 


— 
bo 
eo. 
ow 

— 


Expanding in powers of /, one gets 
a = #,+w,(t—t,)—4ke,(t,—t,)!(t—t,)!—4kc,(t, -t,)+.... (24) 
The change in velocity through the shock may be obtained from a 
combination of equation (21) and the relation 
U,+U, = (1—p?)(w,—we) (25) 
(cf. equation (10)) which holds across the shock. Explicitly we find that 


4(] B*)c, k 


2 - 26 
eral (4) /—t)}'—k i 
——: . 
~ 4k(] 2\e (4 Fis 2 
Be eay5 t, | (27) 


Since |w,+-u.| gives a measure of the shock strength, it follows from 
(27) that, in the absence of a second rarefaction wave behind the shock, 
the strength will decay inversely as the square root of the time. 

Relations similar to those obtained above have been given by Friedrichs 
(2). They are reproduced here for the sake of completeness and to expedite 
the mathematical discussion in Part II where similar results are required. 

For the present problem we require the coordinates (x,,t,) of the inter- 
section of the shock path with the head of the second rarefaction wave, 


namely ia x,+w,(t—t,). (28) 


s 
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A combination of equations (24) and (28) gives the following approximate 


relation for t,: 


tpt, — [Mertech wallet) Ah ev(ty—t)] 
_— fF 4ke,(t,—t,)! 
C(t, —t,.) —(u,—Up)t, —4k?e,(t, —t,) |? i 
ag, (os | (29 
4ke,(t,—t,)* 
since = wf, \ (30) 
Le —Ust, 
and We = Co— Up. (31) 
Using equation (18) and the fact that 
C1—Co p(w, —W,) ee = (Uy— Up) (32) 
a > l—p* _ 


(cf. equations (10) and (11)), the above relation for (t,—t,) may be written 


—¢. — |’ at tolu*t,—t.)—4h(t,—1,)]? 
: ; 4k(t, —t,)! 


r 
r 


(33) 


Interaction of the shock wave with both rarefaction waves 

Beyond the point (x,,/,) of the shock path there is no longer a region of 
uniform state behind the shock. To the present approximation, however, 
the two simple waves are not influenced by the shock. Accordingly the 
shock velocity is given by equation (13), and the flow on either side is that 
arising from the two rarefaction waves defined by 

rL—2, w,(t—t,), (34) 

X—X, w,(t—t,), (35) 
which may, as before, be regarded as equations for the shock path in which 
w, and w, are functions of the time. 


Subtracting equations (34) and (35) we get 


L,—2X,y (w,— w,)(t—t,) + w,(t,—t,), (36) 
and a combination of equations (13) and (34) gives 
| l 
{w,(t—t,)} 3(w,+w,) 4 — (w,—w,)* 
dt 8c, 
lw, l bi 
or (¢ a “1 4(w,—w,) 4 (w.—w,)?. (37) 
dt ailing Sc, 


Finally, substituting for w, from equation (36), we derive the following 
equation for w,—w,: 


d (t—t, 
Gi 


(t (w,—w,)} = 3(w.—w,) + (w, (38) 
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the solution of which is 


W,—W, et arenes » (39) 
2t—t,—t,+ K(t—t,)*(t—t,)! 
; t,—t, (t.—t,|4 
re K s__t — 2{=—_—“"}*, 40 
where k(t,—t,)(t,—t,)! \t,—t,] (40) 
and is determined by the condition that, at t = t,, w, = w, and w, takes 


the value given by equation (21). 
An expression for w, can now be obtained from a combination of relations 
36) and (39); thus 
" Acy(0 t.) 4 Mehr (41) 
. 3-—4-44-Be—the—aF * 4— 
Finally, the equation of the shock path is given by equation (34) using the 


above expression for w,. Explicitly this is 


4c, (t—t,)(t—t, t,— tz, 
y x. l 's) r) | v. vy (¢ -t,) (42) 
2t—t,—t,+ K(t—t,)*(t—t,)* t.—t,. 
2 t—t, 27 + K[ T?—}(t,—t,)?}! 
where Fr t—43(t,+-t,). (44) 
Asymptotic behaviour 
Equations (10) and (39) show that, for large values of t, 
4c,(1—p2 t 
ie A Cal Bw), r) (45) 


(2+ K)t 
According to the present analysis, therefore, the strength of the shock wave 
decays inversely as the time, or equivalently, inversely as the distance from 
the origin. The coefficient of ¢~! in (45) is strictly an approximation to the 
coefficient in the true asymptotic expansion of (u,—u,). Subject to a similar 
qualification, an asymptotic representation of the shock path is found to be 


Irom (43) 


oa Sethe, (ee, Se he Se (46) 
5 1—t,  2+k 27 (2+ Kk)? 


Thus, whereas for a shock wave interacting with one simple wave the shock 
path approximates ultimately to a parabola (see equation (24)), equation 
16) shows that if there is a simple wave on each side of the shock, the shock 
path is approximately hyperbolic at large distances from the origin. 

Of particular interest is the coefficient of 7 in relation (46), namely 
4c, X.—xXx 


- ~, 47 
24+K  t-t -— 


8 r 


1 


since this is the asymptotic velocity of the shock. Correct to the first order, 


an explicit expression for U,, in terms of a is easily obtained as follows. 
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From equations (17), (22), and (33) one finds, after some reduction, that 
—t 


_ &+e 





r ' : 
aE = ae +10), 48 
t,—t,. 7 7 . 9 9 
= tg (+ 8") + 32(1— 2’) + O(0), (49 
\t,.| 2a «(8 7 
—I 
where T= “i hs (50 


If these relations are substituted in equation (40) the following expression 
for K is obtained: 
K = 2—4y?o(1—7-)+ O(o0?). (51 
Equations (47) and (51) now give, correct to the first order, 
oe 4c, ae t,—Usty (52 
2+-K t—~ 


or, substituting for o and 7 from (18) and (50), 








r 


U,. = ht ee. (53) 
(1 —p*)(t,—t,) 
Finally, in terms of the initial velocity of the shock given by (14), 
. r , Uy—U,g t,+1 
UU, = U4+3 4 2 "4+ O(a, —a,)". 54) 
" 2(1—p?) t—t. i( 1 2)*} (5 


As a simple check on this result, we may consider the two limiting cases 
t, > 0, t,—> 0. Equation (54) then gives, to the first order, 
. , , Uy—U; = 
om U4. .4., (55) 
3/1 —*) 
2(1—p 
which, by virtue of relations (14) and (32), implies that 
U, =Cy—Uy OF C,—U, (56) 
respectively, as one would expect. 
Since the shock strength tends to zero at infinity its velocity tends to 
become sonic relative to the motion of the gas. Consequently we must have 
Un = 6.4%, = We. (57) 
If now the value of U,, given by equation (53) is substituted for w, in 
equation (6), we get 
- (6), we g U,t.—Ugt, 


8 (58) 


aes — 


t,—t 


Ss r 
Thus, to the first order, the particle velocity in the neighbourhood of the 
shock tends to the mean value of the piston velocity during the interval 
between the initiation of the two rarefaction waves. 
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Using equations (7), (32), (53), and (57) a similar result may be obtained 
for the sonic velocity, namely, 


Cc ‘2"s Es (59) 


Conclusions 

The preceding analysis has been carried out on the assumption that the 
first interaction occurs between the shock and the simple wave ahead of it. 
However, the basic formulae remain valid if the suffixes r and s are inter- 
changed (implying also an interchange of suffixes 1 and 2) so that the state r 
refers to the conditions behind the shock. The final results therefore also 
remain valid in the general case if we agree to the convention that suffixes 
yand 1 refer to the conditions on that side of the shock including the simple 
wave which first interacts with it. Since an interchange of suffixes leaves 
equation (54) unchanged, this expression for the asymptotic velocity re- 
mains true in the general case without modification. Consequently there 
will be a net increase or decrease in the shock velocity according as the time 
interval between the initiation of the shock and the second simple wave is 
greater or less than that between the initiation of the first simple wave and 
the shock. For the critical case, ¢, t., the present theory says only that 
the net change is at least of the second order. In this case, equations (48), 
49), and (50) show that ¢, and ¢, differ only by a term which is O(1). 
Equation (24) then shows that the change in shock velocity between (2, t,) 
ind (x5, f) is only of the second order. It can also be deduced from equation 
12) that, for ¢ > t,, the velocity is a monotonic function of t. Thus at no 
time during the progress of the shock wave does its velocity differ from its 
original value except by a small quantity of the second order. 

Further conclusions may be obtained from (54) concerning the possibility 
of the interaction of successive shocks formed by repeated impulsive varia- 
tion of the piston velocity. The most interesting case occurs when the 
velocity of the piston repeatedly alternates between —w, and —w., so that 
a succession of shock waves is formed, each with the same velocity and 
separated by simple waves. Thus, in Fig. 3, the two shock waves will 


eventually converge and interact, or diverge, according as 


To determine whether the shock waves ultimately converge or diverge 
in the critical case, 


requires a more accurate analysis. This is the situation envisaged in Fig. 2, 


092,20 
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in which the piston motion is an impulsive oscillation, the velocity changing 
repeatedly from —u to u. In this particular case we may say simply that, 
throughout the motion of the shock waves, their velocity is cy, the velocity 
of sound in the undisturbed gas, correct to the first order. Ultimately their 
strength decays inversely as the distance from the piston. 











Part II. Two-DIMENSIONAL STEADY FLOW 


An approximate procedure similar to that used for one-dimensional flow 
may also be applied to steady supersonic flow along an indented wall, the 
indentations producing rarefaction waves and shock waves alternately. As 
before, it is sufficient to consider a single shock wave bounded on either side 
by a rarefaction wave. The origin of coordinates is taken at the angle of the 
wall initiating the shock, the two rarefaction waves are centred at the points 
(x,,Y,), (%,,Y,) and first encounter the shock at the points (2,, y,), (2.4) 
respectively (Fig. 4). 

Let the suffixes r and s refer to the upstream and downstream sides of the 
shock respectively, let g be the magnitude of the flow velocity, and 0, f the 
directions of the flow and the shock line measured relative to the positive 
x-axis. Then if 

M-1 = sina (1) 
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ging the following relations hold (2): 
that q, cos(B—6,)  p?sin*(B—6,)+-(1—p?)sin@a, 
oa gq,  cos(B—8@,) sin(B—8@,)sin(B—8,) 
their ; ? J 
sin(B—6@,)sin(B—6@,) (2) 
p?sin?(B—6,)-+ (1—p?)sin2a, ™ 
y } Shock path 
Fic. 4. 
From these the following expansions may be derived: 
(24? 1 +-4?)(#2+-p?)t, P , 
ae M6. 0.) ( r)I ro B Mt, 0.—OP+..., (3) 
1—p* (1—p*)* 
| flow 1-2) ° 
oe 1 agy2 
|. the tan(a«.+-6.)—tan(a 0.) uw Ww, n > r) (6, 6,) 
2 
As f 
+ aiilie (1+ #2)(1+-w?) ‘ . 
r sid "ft (f+ p?)+w,(1+4)}(6,—86,)?+ (4) 
of the (1—p*)° 
oints , : l ° 9 = 
tan B Liw tw.) =~ [1—#—2t, w,|(w,—w,)?+.... (5) 
Pos Yo A 8t.(1+-w) 
where f tana (6) 
»f the * 
and UW tan(a+@). (7) 
b the : 
sitive \sin Part I, to within terms of the third order, we may neglect the entropy 
change across the shock front and the influence of the shock wave on the 
(] ipstream flow. To the same degree of accuracy equations (3) and (4) also 


) hold in passing through a simple wave from a state r to a state s (1). 
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Since we are concerned only with a second-order theory, the coefficient 
of (w,—w,)? in (5) may be replaced by its value in the constant flow region 
immediately preceding the shock in which 0, = 0,, «, = a, (say). Thus 


, l ; 
l 4(w,+w,) + — (w,—w,)?, (8) 
8a 
t,(1--w? 
where a= — i v) (9) 


1—t?—2t, w, 
Equation (8), together with the equations of the characteristics in the two 
simple waves, namely 
y—Y, = w,(x—z2,) | (10) 
y¥—y, = w,(x—z,) }’ 
now replace equations (13), (34), and (35) of Part I and the previous results 
may be taken over directly provided that we change the symbols 


t, #, wv, U, t, 


of Part I to the symbols 2, y, w, U, a, 
respectively, of the present theory. 

In particular, the upstream rarefaction wave first meets the shock at the 
point (x,, y,), the coordinates being given by the intersection of 


y = Ux ,+4(w,—w,)-4 o> —W,)? 12 11) 
y x w,+3(w.—w,)- aq (ts v1) ( 
and y—Yy, = w,(x—z,). (12) 
Thus Ly —— 8| Yr a wy 2,| - 82, | tan 9, na wy | - (13) 
4(w.—w,)- (w.—w))*/a do(4-+-oc) 
where o = (wy—w,)/a. (14) 


Interaction with the first rarefaction wave 
With the appropriate change of symbols equation (23) of Part I gives, for 
the equation of the shock path between (x,, y,) and (2., ys). 


4ka(x,—x,)*(a—a,)! 


n y+-w,(%—x,) ———— ron (15) 

y Yr o( ) l kf (2, 2,.) (x—a,)\# 
Y,+w,(x—2x,)—4ka(x,—2,)(a—a,)!—4k?a(x,—--a,)+..., (16) 
where elena, (17) 


The coordinates (,, y,) of the intersection of the shock line with the head 
of the second rarefaction wave are obtained from a combination of (16) and 


(18) 








Thus 


This 


Inte 
Th 


equa’ 


wher 


and 
(cf. € 
Fi 
theo 
tion: 
quar 
weig 
rare’ 


A 
is pe 
ang! 

lr 
of tl 
stea 


T 


and 


cient 
2101 


1us 


> tW 


(10 


sults 


t the 


(1] 


(12 


q fi yT 


(16 


(17 


lead 


and 


(18 











DECAY OF SHOCK WAVES 


Thus we find that 


(19) 


Y,—Y. W(X, —2Z,)— 4kh*a(x, -- a] 2 
a 
4ka(x,—x,)* 
This relation corresponds to equation (29) of Part I. 
Interaction with both rarefaction waves 


The behaviour of the shock beyond the point (x, y.) may be obtained from 
equation (43) of Part I, namely, 





y- ) 8. Iie... 8 
y Y, Y. | Y. Y, T | _4a T ' 4%, x,) on (20) 
2 L,—2 2T + K|[ T?—}(x,—z,)?}! 
where K Tet ig (tard (21) 
h(a, L,)*(Xy—2%,)* \a— x,J 
and 7 o— 4(x, T Z,) (22) 


(cf. equations (40) and (44) of Part I). 

Further properties of the flow are similarly deduced from the unsteady 
theory. Thus, for example, the strength of the shock is inversely propor- 
tional to the distance from the origin and, as can be shown directly, all flow 
quantities in the neighbourhood of the shock tend asymptotically to the 
weighted mean of their values in the constant flow regions bounded by the 
rarefaction waves and the shock. In particular, to the first order, 


| Xo -F 4.) — x,( Xy +-8,) 


B. ¥. +6, et —, (23) 
X,—2%, 
—2x 
4 Oo — 2, 8 (24) 
oe — Ly 
XL.Jo—2X,G " 
q 8s 12° rQ (25) 
x,—i, 


NOTE ADDED 10 DEC. 1951 


A case of practical interest is that of steady flow past a rough wall whose mean shape 
is parallel to the mean velocity and which has indentations of equal length all making 
angles +4 (say) with this direction. 

In this case the asymptotic values of the flow quantities B, 0, q in the neighbourhood 
of the shock are, to the first order, the values of a, 0, q in the equivalent problem of 
steady flow past a smooth wall. 


The coordinates (z.. y,) now satisfy 


U, Ys J P 
ind 6, 6, = ¢. 


It is convenient to regard the flow upstream of the point (x,,y,) as steady and 
} re JT . 
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parallel to the x-axis. If the suffix 0 now refers to conditions in this region, the detailed 
flow is described in the following formulae: 


2(1+8)*¢ 





ty( 1+?) 
a ; ail 
1 — 37 (9 bis 
Wy—W, 
one a: (14) 
a 
Sx, t , of 1—p? | 
‘ 4+ “3la gt1y]- 3 bis) 
weal: 2\1+#)2 1; (13 bis 
The qua ion of the shock path between (a iY) and (95 Yo) is 
tha(x, } x) E(x 1 a)2 
y Y. W(X ee Ee i i” ’ 5 biel 
. l ky (a, +x.) /(a- x,)}4 (15 bis 
o 
where k P 
tt+o (17 
(1+)? to( 1 +42)? 
and Ws ‘ d+ oy re 21162, 
. i—" (1—p*)* 


At the intersection of the shock with the head of the second rarefaction wave 


2w.x,— 4k?a(x,+2,)]? $ 
ils tha(x, x,)3 : (19 bie) 
Finally, beyond the point (2, ¥2), the equation of the shock is 
4a(x? — a2) : 
| Wet Se Kaa mn 
. 2x, (ary+-a,\4 ; 
where K hie, ae, a) "\s. a (21 bis 
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ON THE EQUATION OF DIFFUSION IN THE 
ATMOSPHERE 


By E. KNIGHTING (Meteorological Office, Shoeburyness) 
[Received 11 September 1951] 
SUMMARY 
An analysis is made of the diffusion equation in two dimensions with the diffusion 
efficient proportional to a power of the height above a given datum level, nomi- 
nally the surface. This equation is used in meteorological problems concerning the 
wer atmosphere ; solutions are given under the usual boundary conditions met with 
these problems and solutions of particular problems are indicated. Special atten- 
tion is paid to the possible values of the index in the power-law diffusion coefficient, 


and conclusions drawn regarding the applicability of the solutions in various cases. 


1. Introduction 
Ix considering turbulent transfer in the vertical in the atmosphere the 
differential equation often used is 


Cc ov dv ov ov 
az'-* ; —+u—, 
C2 dt ot Ox 


Cn 


where x, z, and ¢’ are the horizontal, vertical, and time coordinates, u the 
horizontal velocity, and v is the entity under discussion; a and m are con- 
stants in a particular meteorological situation. Sutton (1), inter alia, has 
discussed this equation in the case when év/ét’ = 0 and uo 2”, following 
the lines of the classical treatment of the conduction of heat. The equation 
in which « = 0 has been treated in special cases by various writers, e.g. 
Frost (2), Knighting (3). Kohler (4) has discussed the equations for 
0<m 1 in a very thorough treatment. 

It has been conjectured that near midday ona clear summer day the index 
in the diffusion coefficient for heat lies between about 1-5 and 2 (Sutton (5)). 
Experimental work carried out at Rye (Sussex), and described by Sir Nelson 
Johnson (6), provided data concerning both temperature and humidity dis- 
tribution with height; these data were examined in order to determine what 
relation, if any, there is between the diffusion coefficients of heat and water 
vapour in certain meteorological conditions. During the course of this work 
difficulties were encountered concerning the possible range of values that 
the index in the power law may take. 

We give here an account of the solutions of the equation with the boundary 
conditions usually met with in meteorological problems and indicate how 
specific problems may be attacked. A generalization of the equation is made 
in order that the two cases quoted above may be dealt with at the same time 
ind the possible values that m may take are stressed. 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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KNIGHTING 
2. Some mathematical formulae 


For convenience some of the mathematical formulae are collected here 
for reference. The p-times Laplace transform is used, defined by 


fp) =p | exp(— pt) f(t) dt, (1) 
0 
and we write f(t) >f(p). (2 
Two particular results required are 
f(t—h) > exp(—hp)f(p) (t > h), (3 
L. x? + B? x8 
5i°8P( > )4n(5) > pl, («p*)K,,(Bpt) (n > 0), (4) 


where /,,(x), K,,(a) are the modified Bessel functions of order n. In particular, 


nO n+2j 
= 1% rj (n+j+1)’ 9) 
) 


x 


. P - 
a" (2) = 2° | wenn (6) 


0 


This equation holds for x = 0, provided that the integral converges. Also 


K,,(x) = K. n(X), (7) 
and . 292 (2) = 2°f, _,(2), . x"K, (x) = —2x"K,,_,(2). (8) 
dx dx ” 


3. The equation and the boundary conditions 
The equation to be discussed here is 


¢ ov ov 
—{az'-* —) — kz —, (9) 
02 Cz ot; 


where ¢, is either 2 or t’, k and / are constants. Clearly, with proper choice 
of k and /, equation (9) covers both cases mentioned in section 1; moreover, 
if the conjugate power law (J = m) does not hold but w is still proportional 
to a power of z, then this equation provides the solution for steady motion. 
The substitution 

at, 


t = _* at’ in the time equation, 


ax. R 
—-— in the space equation, 


k 


removes the parameters and (9) becomes 


3 gl—m ov eee ov (10) 
Oz OZ et 
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ON THE EQUATION OF 


In meteorological applications we are usually given the distribution of v 
with z at f 0, and the variation of v with ¢ at the surface z = 0. In 
addition, the value of v must remain finite everywhere, and in particular, 
as z> 00. Other boundary conditions which may arise, such as a given 
value of v over a finite strip, as in evaporation problems, may be equally 
effectively dealt with by the method of solution adopted here. In some 
problems the boundary conditions may be replaced by specifying v at a 
given height h; Jaeger (7) has dealt with problems of this type. It is sup- 
posed that the value of v at great heights is little altered by the turbulent 
motion and in the solution v will be measured from its constant value at 


z+ 0. The boundary conditions are therefore 


v>f(t) asz>0 (11.1) 
vy > dz) as t> 0 (11.2) 
v—> UV as 2—> oO. (11.3) 


The solution breaks up naturally into two parts. If U satisfies equations 
(10) and (11) with ¢(z) = 0, and V satisfies the equations with f(t) = 0, 
then v = U-+V satisfies the equations as given, as is easily seen. 

In the solution of these equations the following abbreviations are used: 


m 


n , 12.1 
l+m-+1 ( 
2r = l+m-+1, (12.2) 
9oi(l+m+)) yh ah yt 
2 me aagph we Scene sale ale, (12.3) 
l+-m-+1 r 
Qyil+m+Dpi y"p* 
Bp! : - = ——_., 12.4 
e f l+m+1 r ( 


where y is a value of z to be defined later. 


4. Solution of the equations for U 
Applying the transformation of equation (1) to equations (10) and 
(11) with ¢(z) = 0, we find the system of equations 


(2! mol } p2U, (13.1) 
C2 Cz | 

U+f(p) asz>0, (13.2) 
U > 0 as z-> 0. (13.3) 


Jaeger (7) has developed a solution of these equations with | = m, 
f(t) = a constant, proceeding in a rather different manner from that given 
here. The basic solutions of equation (13.1) are 


U = d"I,(A), U = A"K,,(A). (14) 
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Equations (13.2) and (13.3) show that only the K-form is applicable and 
that n must be positive, so that 


U = Ad*K,,(A), (15) 
while from equations (6) and (13.2) 
F(p) = A2"T(n). (16) 
Equations (6), (15), and (16) give 
7p 1 ¢ 2*"p . 
U = = »xp| —p— ate ip. 7 
Fin) | exp/ p on) J (p) dp (17) 
0 


We may apply equation (3) directly to this expression, giving 


x 





l i. ~2r 
Uv = n—lex , | eee lp. g 
rin) p” Texp( as( an) dp (18) 
2"/4rt 


In equation (15), A” is effectively z!” and to ensure that U is finite as z> 0 
it is necessary that m > 0 and we have 


m ; 
nies ae ee (19.1) 
l+-m-+1 
m > 0, i--m-i 2, (19.2) 


Normally, / is required to be positive or zero, otherwise in the steady state 
case the wind speed at the surface becomes infinite. 
The flux at height z and time ¢ is given by 


ov ’ 7 
flux oe - oF (z, t) o F(z), (20) 

Oz 
where cis to be determined in accordance with the meteorological variable 
denoted by v (if v is the temperature, o = pc, where p is the air density 
and c,, the specific heat at constant pressure) and we write F(z) for F(z,t) 
to avoid continually writing F(z, p) in the sequel. Considering the flux of U, 
we have, by equations (15) and (16), 

3 eu ri m . 


F(z) ~l—m “s A! a _* n(A)p"t(p). (21) 


Cz C(n) n(2r)?” 


using equation (8). Again, since n is usually fractional then, for real values 


of the flux, r must be positive with m and n. The flux at z = 0 is given by 
3 r(1 z 
F(0) m U ”) (or) 2nnnf(p), (22) 
nT(n) 


' m(2r)-" d [ , 
F(0) t—r)-"f(r) dr. (23) 
n(n) dt } \ mi ye 
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By writing equation (23) in this form some difficulties are avoided at the 
upper limit of integration arising from the presence of a term (t—r)-"-! in 
the integrand. 

A special well-known case in the conduction of heat arises by putting 
|= 0,m = 1,sothatn = sandr = 1. Setting f(t) = constant = 7’, then 
equations (18) and (23) reduce to 


U T erfe(z/2t!), 
F(0) =, 
(art)? 
which are well-known formulae (Carslaw and Jaeger (8)). 

In a practical case it is better to use equation (22), expanding f(p) in 
negative powers of p and interpreting by the rule 


t” 


['\(n+1) sda 


= 9 


For example, if f(t) = sin wt, f(p) = wp/(p?+w?), 


F(0) ml (J M) (2p) npn fy = F o 
nI'(n) | p* 

F(0) mt M) y(2r) on{ Ra n al n , 7} 
nI'(n) ir(2—-n) T(4—n) | 


The most important equations are (18) and (23), with the condition that 
nis positive. If m exceeds unity it appears from equation (10) that év/éz is 
infinite at the surface, z = 0, which may or may not be acceptable from a 
meteorological point of view, and indeed is a point which cannot be tested 
experimentally. The flux at the surface remains finite in this case. 


After a sufficiently long lapse of time equation (18) may be replaced by 


: E> ¥ ual 
l ; p” exp( pdf -—~ } dp. (24) 
I'(n) P 4r“p 


0 


But a warning is necessary. The greater part of the contribution to this 
ntegral comes from values of p near the lower limit, owing to the presence 
(the exponential term, and it is necessary to show that the neglected part, 
the integral from 0 to z2"/4r2t, is ver v much smaller than the infinite integral. 
For diurnal heating in the period from dawn to dusk an approximation to 
the form of f(t) is sin wt, where the period of the wave is 24 hours. Between 
dawn and midday only one-quarter of a wave has passed and the 


asymptotic form of equation (24) may not be a good approximation. 
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5. Solution of the equations for V 


Applying the transformation of equation (1) we find the system 


C oV oa 
= (2-* us ) = 2'p{V —d(z)}, (25.1) 
ez oz 
V>0O asz> 0, (25.2) 
V+0 asz>o. (25.3) 


To solve, we write the solution with ¢(z) = 0 as w(z), where 
w(z) = w(1) = Ad", (A) (0<z<y,ie.0<A< yp) ( 

= w(2) = BAK, (A) (y¥<z2<0,u<A< oo), (26,2) 
Be 


for these two solutions satisfy the boundary conditions, equations (25.2) and 
(25.3). In addition we require that 
w(l)=w(2) atz=yorA=uz, (27.1) 
’ Agn(D 
i =e" ats= 9, A= p. (27.2) 
az 2 


The last two equations serve to define A and B in (26.1), (26.2); it is easily 


found that ptr K, (w)D,,(A) dy 





(1) = 
w(1) yim des 
9 pl j mr” I, (u)K,,(A) dy 
w(2) = x ——" _" , 
y m du 


and it follows that the solution of equations (25) is 


oo 





V(z) = - | pysty) ae | | py’b(y)w(2) dy. (28) 


0 2 


The first integral, Z(1) say, reduces by use of equations (12) to 


0% 


L(1) = = | a"B'-"b(y) pI, (Bp*)K,, («p*) ap. 


0 


The second integral is similar with J, and K,, reversed, the limits being (a, 0). 
Using equation (4) we find 


a sh | exp(— B ae nefs( (B)I, “() dp, (29) 
where ¢(y) = (8). : 


The proviso in the interpretation of L(1) is again that n must be positive. 
The flux is most easily obtained from equation (29) and gives 


F(0) = — —m 


Tint 77 | bald . | By(Byexp( —") dp, (30) 


and we see that r must be positive, and hence m must be positive. 
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Turning again to the special values of the parameters for the conduction 
of heat we recover the familiar formula 
: Yi 
F(0) ——... 
(zt)? 
It is interesting to calculate the two fluxes when f(t) = ¢(z) = 7’, for then 
is everywhere constant and equal to 7’, and there should be no total flux. 
From equation (23), with z 0, we have 
; om(2r)-2"t-" 
flux for l ) eaeoaiarti 
n(n) 
while, from equation (30), 


flux for | __om(2r) 2n¢ n 
| n(n) 


showing that the two fluxes are equal and opposite. 

For large values of t, V tends to zero whatever ¢(z) may be, provided that 
it is such that the integral in equation (30) exists, which implies that after 
, long time the distribution of v at ¢ = 0 is immaterial, the transfer being 
determined by the impressed variation at z = 0, andv ~ U. 

The solution for V is very well adapted to calculation if ¢(z) = ¢(a) isa 
polynomial, for then we may expand J, (af/2¢) in powers of 8 by means of 
equation (5), and perform the successive integrations. If z is small, a first 
approximation is given by the first term of the Bessel expansion, and 


writing . 
ing A( 2) ub( x) - | x") 


x 


then V ore | exp/ Ai) Mee) age 


(4¢)"*11D(m-+-1) | 4 
0 
If 6(B?) b a, B®, 
then fa > — . U Sis. e"exp| — af 
La (4t)"-§ P'(n+-1) 4t 


and, in particular, if d(z) = z = (a?)"/" then 


y?\(m—I)+m+D) PY] 4-1 /(l+-m-+1)} x” 
(*) | + oun — =) 
4t 1+-m/(l4+-m-+1)} 4t 
This result has been given by Frost (2) for 1 = m, except that Frost omits 
the term exp(—.«2/4t) which is nearly unity for small values of z. 
6. Discussion of the solutions 
The solution of equations (10) and (11) is now complete and 


v= U+V. 
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Only the formal solutions have been dealt with and where infinite integrals 
occur it has been assumed that they exist; e.g. in equation (18) it is assumed 
that f(t) is such a function that the right-hand side of the equation exists, 
The proof that equations (18) and (29) do satisfy the differential equations 
and the boundary conditions can easily be supplied. No proof that these 
solutions are unique is attempted and indeed such proofs would be out of 
place in a practical paper of this kind. 

Some remarks upon the applicability of the solutions are not out of place, 
It should be recognized that the diffusion law has been assumed to hold to 
great heights, indeed to z—> oo. In practice it is scarcely likely that this is 
true and solutions are approximately true in a layer extending up to some 
height h; it should also be noted that the larger z is in the interval (0, h) the 
poorer the approximation, and great caution is necessary in applying asymp- 
totic forms in z, lest the asymptotic form be not a good representation of the 
facts in the upper part of the interval (0,/). The caution necessary in using 
asymptotic forms in ¢ has already been noted. 

A second important point is that if the experimental results are analysed 
by means of the solutions given above and lead to values of m which are 
negative, then the solutions are not applicable. Alternatively, if experiments 
result in values of m that are negative, no use having been made of the 
solution, then the solutions are invalid. In either case the diffusion equation 
cannot represent the motion. It is concluded that, with the boundary con- 
ditions as stated, the index in the power law representing the diffusion 
coefficient is less than unity, if the differential equation is to represent 
adequately what happens. 

For large values of z a solution for U can be obtained from the asymptotic 
expansion of the Bessel function in equation (15), and it is found that there 
is no restriction on the values of m. This asymptotic form has been used in 
the analysis of temperature changes in the lower atmosphere, assuming 4 
sine wave of temperature at the surface; the latter assumption is only valid 
on days of clear weather. The phase shift of the time of occurrence of the 
maximum temperature at a given height can be calculated, and the observa- 
tions show that the change of phase is slow, roughly as 2°! or z°*. Ifm > 0 
the phase shift cannot be less than 2°, so it seems that the power-law 
diffusion coefficient K = az!-™, which is a great improvement on K = con- 
stant, for the transfer of matter and momentum, does not apply equally well 
to heat transfer in the atmosphere. This may be because the transfer of heat 
involves both free and forced convection, and fine weather is favourable for 
free convection whereas the transfer equation that we have used is really 
only applicable to forced convection; we must not make the observations 
for the case of free and forced convection together fit into the solution for 
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forced convection alone. The fact that the asymptotic solution does agree 
with observation when m < 0 but that the solution fails near the ground 
unless m > 0 may indicate that the index m varies with height. I am 
ndebted to Professor O. G. Sutton for drawing my attention to the points 
raised in this paragraph and for his general criticism. 
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ON THE INTEGRATION OF SOME VECTOR 
DIFFERENTIAL EQUATIONS. I 
By G. N. WARD (Department of Mathematics, The University, Manchester 


[Received 29 November 1951] 


SUMMARY 

By using a suitable vector identity, which is derived in section 2, it is shown that 
the system of inhomogeneous first-order vector differential equations (1) given below 
can be solved directly in vector form for suitable boundary conditions. Such systems 
of equations occur in a number of branches of mathematical physics. 

Two main cases arise, depending on whether the system of equations is elliptic or 
hyperbolic. The solution in the elliptic case is straightforward, but convergence 
difficulties occur in the hyperbolic case which are overcome by using Hadamard’s 
finite part technique for dealing with fractionally infinite integrals. 

The results are given in invariant vector forms that do not depend upon any 
particular choice of coordinate axes. 


Notation 

GiBBs’s dyadic hotation is used for the representation of a second-order 
tensor, and ¥.f denotes the inner product of the tensor Y and the vectorf 
(cf. C. E. Weatherburn, Advanced Vector Analysis (Bell, 1928), chap. v): 
in cartesian tensor notation, if f, 8, and WY have components /,, g;, and ¥; 
respectively (7, 7 = 1, 2, 3), then the last equation of the system (1) below 
can be written as . 


3 Yuh 


1 


9; 


In the case under consideration, all the components a oF are constants, and 


A ; A 


ij ji* 


1. The vector equations 
In a number of branches of mathematical physics, systems of vector 
differential equations are encountered that can be written in the form 
¥AE = P, V.¢ = Q, ¢= ¥.f, (1) 
where P and Q are respectively given vector and scalar functions of position 
in three-dimensional space, Y is a constant symmetrical tensor of the second 
order, V is the gradient operator, and f and g are unknown vector functions 
which have to be determined from boundary data given on some surface. 
For example, the electrostatic and magnetostatic equations for homo- 
geneous anisotropic media are respectively (1) 
VAE a ¥.D = g, D = «.E, 
and VAH = J, ¥.B = p*, B= «.H, 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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where, in suitable units, E, D, H, B, J, and p have their conventional 
neanings, J* and p* are respectively fictitious magnetic current and charge 
lensities, and € and pw are dielectric and magnetic permeability tensors 
1. p. 11). 
Again, the linearized equations of motion for the steady flow of an 
viscid compressible fluid can be written as (1) 


VAN a Vw Q, Ww (I— UU /a?).v, 


t 


where V is the perturbation particle velocity vector, ¢ is the vorticity, Q is 
the source density of volume, U is the constant velocity of the undisturbed 
stream, a@ is the acoustic velocity, and I is the unit tensor or idemfactor. 

In the special case when P = 0 in the system (1), f is the gradient of a 
scalar potential function, ¢, and the system can be reduced to a second-order 
partial differential equation ford. This equation can be reduced to Poisson’s 
equation or the inhomogeneous two-dimensional wave equation by suitable 

hoice of rectangular cartesian coordinates and an affine transformation, 
n which case the solution can be obtained by standard methods. 

For another special case, when Q = 0, and g = constant xf, a vector 
solution has been constructed by Stratton, using a vector analogue of 
Green's theorem (1, p. 250). 

It is shown below that vector solutions of the system (1) can be con- 
structed in the general case by using a vector identity which was derived 
by F. Ursell and the author for the purpose of proving uniqueness and 
reciprocal theorems (2). Extreme mathematical rigour is not considered to 
be the most important part of the investigation and is consequently omitted. 
The functions under consideration are assumed to be sufficiently regular 
everywhere to ensure the validity of the transformations used. In any 
practical applications of the results, it is necessary to check this point, but 
conditions are so varied that it is better to do this as they arise rather than 


to try to lay down sufficient conditions here. 


2. The integral identity 
The vector identity required for the construction of solutions of (1) can 


be derived by applying the divergence theorem to the tensor 
éf—if.gsl; 


if nis the outward unit normal to the surface S enclosing the volume JV, 


this gives 


{ (fg.n—1f.gn) dS = [ ff£V.g—8A(VAf)} dV. (2) 
S V 


5092.20 _£ 
F f 
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Now put f = f,+f,, etc., and f = f,—f,, etc., successively in (2) and 
subtract, which gives 
[ (f,8..n-+f,6,.n—f,.g,n) dS 
- 
7 | if, V-8o+f, V.81:—81 A (VA f.)—82 A (V Af,)} dV, (3 
Vv 
since f,.g, = f,.¥.f, = ¢,.f,. This is the required identity. 


3. The homogeneous equations 

For part of the following analysis it is found convenient to introduce 
rectangular cartesian coordinates x, y, z, with axes such that the tensor ¥ 
has only diagonal components. Ifi,j, k are unit vectors parallel to the axes, 
Y can then be written in dyadic form as 

Y = a, ii+a,jj+a, kk. (4) 

The character of the solution of the system (1) depends upon the signs 
of the coefficients a,, a,, a, in (4). This can be seen by considering the 
homogeneous equations 


VAf = 0, V.¢ = 0, ¢= ¥.f. (5) 
The first equation of (5) is satisfied by writing 
f <e V4, (6) 


and the second and third then give 
V.(%. Vd) = 0, (7) 

or, in the special coordinate system, 
Ay Pyzt+Ae byy +43 $2. = 9, (8) 


where suffixes denote partial differentiations. 

Equation (8) shows that if a,, a,, a, have the same signs, the system (1) is 
elliptic; if not, the system (1) is hyperbolic. These two cases have to be 
treated separately. 

It is assumed below that none of the a’s is zero or infinite. 


4. The elliptic case 
4.1. Consider first the elliptic case, and take 
1, Az, Az > 0. (9) 
No loss of generality is involved in this assumption, for if the a’s are all 
negative, a change of sign in Y, and f and P, say, restores the condition (9). 
A fundamental singular solution of (8) is then 


_ {(z—a,)? | (y—m)? | (@—%)\-4 (10) 
eed ame ‘ “i ii | : 


— —}— — 
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If the position vectors of the points x, y, z and 2, y,, 2, are respectively 
R and R,, and ¥~! is the reciprocal of ¥ (so that ¥ Y-! —Y-1 Y — J), 
then (10) can be written 
¢, = {(R—R,).¥-!.(R—R,)}* = Rz', say, (11) 
which gives the fundamental solution in a form independent of the special 
choice of axes: R, can be called the ‘elliptic distance’ between the points 
Rand R,. The corresponding values of f and g are 
f Y-! (R—R,)/R3, ¢, = —(R—R,)/ RF. (12) 


This fundamental solution (10), (11), or (12) is the analogue of the 


solution of Laplace’s equation (when a, = a, = a, = 1) for an isolated 
pole or sink, and there is a constant flux of the vector g, through a closed 
surface S which surrounds the point R,. If n is the unit outward normal 
to 8, this flux is ‘ ' ‘ 

| ° |} &o.ndS, (13) 
Ss 

and the integral is invariant with respect to changes of S since V.g, = 0, 
R + R,; thus it can be evaluated for the case when S is a spherical surface 
with centre R,. This evaluation is most easily carried out by introducing 
spherical polar coordinates |R—R,|, #, @ with origin at R,, and polar axis 


parallel to the x-axis, say, in which case 


So n dS 








cos’ cos*g _ sin’\]-3 
26 T 2 - 
| |- sin H + sin ? dtde 


As As 


47(a, 4,43). (14) 
But a, 4, a, is the determinant of the elements of Y in the special coordinate 
system, and this determinant is an invariant of YW. Hence, finally, for any 
choice of coordinate axes, 


( $.ndS = —4n(det ¥)}. (15) 


4.2. Now let S be a closed surface enclosing a volume V, let R, be any 
point of V not on S, and let S, be a sphere of small radius with centre R, 
enclosing a volume JV, and lying entirely in V. Then a solution of the 
elliptic system (1) can be constructed by applying the identity (3) for the 
surface S-+S, enclosing the volume V—J,, with f,, 8, as the required 
solution f(R,), &(R,) at the point R,, and f,, g, as the fundamental solution 
f,, 8) of section 4.1. 
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The contribution from the surface integral over S, is 
f | g).ndS + 8| f,Ands, (16 
Si g. 
where the bars denote mean values over S,. The value of the integral in 
the first term of (16) is 47(det Y)!, from (15), remembering that n is noy 
the inward normal to Vj, and the integral in the second term vanishes 
identically, as can be shown by expressing the integrand in terms of spherical 
polar coordinates. Hence, in the limit when the radius of S, tends to zero 
the contribution from S, is 


47(det Y)#£(R,). (17 
The surface integral over S can be written as 
((n. $)fy+ (n/ f) \ Bo} dS, (18 
RS 
and the volume integral is 
| (Qf,—8 AP) dV, (19 
V-V; 


since f, § satisfy (1), and f,, 8, satisfy (5) in V—N. 
Hence, if the volume integral converges in the limit when the radius of 
tends to zero, (3) and (17) to (19) give the identity 
* {nA f(R)} \ (R—R,) 
R? 


17(det Y)!4(R,) | n.g(R)V( >) as 1 dS 4 
v, 
* (R—R,) A P(R) 


Ri dV. (20 





. ix 
| QR)V( 7) aI n 
F 


A corresponding identity for § is obtained by multiplying both sides of 
(20) scalarly by ®. 

Thus, if a regular solution f, g, of the system (1) in V, having specified 
values of n/ f and ‘n.g on S, exists, it must have the form (20). 


r 


In order to verify that such a solution exists, it must be verified that it 
satisfies the system (1) and reduces to the given boundary values on 5. 
That the solution (20) satisfies the system (1) can be verified provided that 
nf, n.g, P, and Q are sufficiently regular, but no verification of the 
boundary conditions is possible, since nf and n.g cannot be specified 
independently on S. 

4.3. That the solution of (1), if it exists, is unique when n / f is specified 
on some part of S, S’ say, and n.¢ is specified on the remainder of S, S” say, 
can be demonstrated as follows. Suppose that two solutions can be found, 
different if possible, and let f, g be their difference; then f, ¢ satisfy the 
homogeneous system (5) and the boundary conditions n A f = 0 on 8’, and 
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n.g=0on 8S”. As in (6), f = Vd, and ¢ must be constant on S’. The 


divergence theorem applied to ¢g gives 


| dg.ndS = | (f.g+4V.8)dV = | (a, ¢2+a,4¢5+a,¢2)dV. (21) 
k i Vv 
Whatever the constant value of ¢ on S’, the left-hand side of (21) vanishes 
since . 
}¢.ndsS | g.ndS = 0; (22) 
S Ss 
ind since the right-hand side of (21) is essentially positive, uniqueness 
follows by the usual arguments. 
Other uniqueness results for more complicated boundary conditions can 


be demonstrated in a similar way. 


4.4. Solutions for the boundary conditions of section 4.3, or for more 
complicated conditions, can be constructed if suitable Green’s functions 
for Sare known. For example, suppose that the solution is required when 
n.g is specified on S, which is the analogue of Neumann’s problem. Then 

vector Green’s function, F(R, R,), must be found, satisfying the system 


and the conditions 


nA\F=0 onS | 
F~f, as |R—R,| > 0 


. (23) 
F is regular in V except at R RJ 


If G(R,R,) is the scalar Green’s function for equation (7), satisfying the 


conditions 


(s () on s&s | 
G~ R>* as |R—R,|->0 , (24) 
G is regular in V except at R R,| 


then F = VG is the required vector Green’s function, and the required 


solution, if it exists, must be given by 


47(det ‘Y)*£(R, ) } n.g(R)F(R,R,) dS 
S 
| Q(R)F(R, R,) dV + | (%.F(R,R,)} A P(R)dV. (25) 
V 
If the integrals exist, a verification that f(R,) satisfies the system (1) and 
reduces to the correct boundary values on S is now possible. 
Vector Green’s functions for other boundary value problems may be 


mstructed from corresponding scalar Green’s functions in a similar way. 
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5. The hyperbolic case 
5.1. In the hyperbolic case take a, > 0, az, ag < 0, and write ¥ in the 
form : 
Y = b, ii—b, jj—b, kk, 
where 6,, b,, 6, > 0. Again no loss of generality is involved by this par. 
ticular choice of signs. Equation (8) then becomes 


by Bus —b, Pyy —bz Pez = 0, (26 
which has the fundamental singular solution 


db {(x Hy)? (yy)? _ (@—%4)? 5 
Po 
\ 


by b, bs 

((R—R,).¥-!.(R—R,)}-? = R,', say, (27a 

for Ri > 0, x < x,, and 
od, = 0 at all other points. (27b 

The corresponding values of f, and g, are 

f, = -¥1.(R—R,)/R3, 8, = —(R—R,)/R} (28 

for Ri > 0,x < x, and 
., = 0, é, = 0 at all other points. (28b 


Besides being singular at R = R,, this solution is singular on the charac: 
teristic half-cone R, = 0, x <2,. The quantity R, has been called the 
‘hyperbolic distance’ between the points R and R,. 

The construction of the solution for the hyperbolic case follows much 
the same general lines as that for the elliptic case, but some complications 
arise from the more singular character of the fundamental solution, which 
causes the integrals in the identity (3) to become divergent. The device to 
be used in overcoming these difficulties was indicated by Hadamard (3), 
who introduced the idea of taking the finite part of fractionally infinite 
integrals. Hadamard’s method has been generalized by Robinson (4) for 
the case of vector functions. The method is now so well known in the 
treatment of hyperbolic equations that no further explanation will be 
given here. 

The flux integral for g,, taken over a surface S which encloses R,, diverges, 
but the finite part can be taken in Hadamard’s sense. It is no longer cleat 
that the finite part of the flux integral is invariant with respect to changes 
of S, but it is not difficult to show that this is so, and S can be taken to 


be a spherical surface with centre R,. The only contribution to the integral 
comes from that part of S for which Rj > 0, 2 <2,. By introducing 
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spherical polar coordinates, as in section 4.1, the flux integral is 


+ r *? [cos?d . 5,/cos’e . sin?o\]-3 . 
¢,.ndS | do | —sin?? Eo. |)! *sind dd (29) 
. " @ | b, b, Is 
. 0 0 
oe cm mT 
‘ COS" sin“< 5 
where V1(9) cot a j + j ) 
- I3 





ind the star in front of the integral sign denotes the operation of taking the 
finite part. This integral can be evaluated (3) to give 

“| $,.n dS = 2n(b,b,b,)t = 27(det ¥)!, (30) 
the latter expression being valid for any system of coordinate axes. 

5.2. In order to construct the solution of the system (1) in the hyperbolic 
case, the identity (3) is applied to the domain V bounded by the charac- 
teristic half-cone R, = 0, x < x,, and those parts of a surface S and a 
spherical surface S, with centre R, which lie inside the half-cone, as illus- 


trated diagrammatically in the figure below. 





As before, take f,, g, to be the required solution of (1) at the point R,, 
ind f,, 8, to be the fundamental solution f,, g, of section 5.1. The surface 
and volume integrals diverge, but the finite parts can be taken (4). 

The contribution from the surface integral over S, is 


f "| g).ndS + 8A “[f ands, (31) 


where the bars denote mean values over S,. As the radius of S, tends to zero, 
27(det Y)'f(R,), (32) 


and the second vanishes ident ically as before. 


the first term yields 


The finite part of the surface integral over the half-cone is zero, and the 
surface integral over S has the same form as that in (18), except that the 


inite part is required. 
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Thus, finally, in the limit when the radius of S, tends to zero, the identity 
(3) gives 
*r inAf(R)} A (R—R 

S) v as _ 


3 
. h 





2n(det ¥)'(R,) = | n.g(R)V Jas . 
J R,, 


(33) 


«0 ‘1\., *¢ (R—R,)AP(R),., 
: | ARV pe) aI -{s = dV, 
J 


* 
and this must be satisfied by any solution of the system (1). 

To show that f(R,), given by (33), is the solution of (1) having the given 
values of nA f and n.gon S, it must be verified that this solution satisfies (1) 
and gives the specified boundary values on S. The former can be done, but 
the latter can be done only if S is everywhere ‘duly-inclined’ in Hadamard’s 
terminology (3, p. 44), or else is ‘characteristic’. In the case considered 
above, S is duly-inclined if the normal n satisfies the condition n.W.n > 0, 
and is characteristic if n.¥.n = 0. The verification of the boundary con- 
ditions is long and difficult, and its possibility here is argued only by analogy 
with Hadamard’s verification in the scalar case (3). When S is duly-inclined 
the problem is analogous to a properly-set Cauchy initial value problem for 
a second-order partial differential equation, and n,f and n.g can be 
specified independently on S; the solution is then unique. 

If S is not duly-inclined everywhere, then verification of the boundary 
values is not possible, and n / f and n.g cannot be specified independently 
on those parts of S which are not duly-inclined. In order to construct the 
solution in this case, suitable functions analogous to Green’s functions for 8 


must be found, and verification of the boundary values may then be possible. 
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ON THE INTEGRATION OF SOME VECTOR 
DIFFERENTIAL EQUATIONS 
II. APPLICATION TO THE LINEARIZED THEORY OF 
STEADY COMPRESSIBLE FLUID FLOW 
By G. N. WARD (Department of Mathematics, The University, Manchester) 
Received 29 November 1951] 


SUMMARY 


It has been shown in Part I (1) that solutions of a special system of vector differen - 


11 equations can be constructed by using a suitable integral vector identity. 
\ special case of the system arises in the linearized theory of the steady flow of an 
iscid fluid, and in this paper the general solutions of Part I are adapted to the 
flow problem. 
Special interest is attached to the particular integral of the equations which gives 
perturbation velocity in terms of the source and vorticity distributions. These 


easily modified to give expressions for the velocity fields due to source layers, 


tex sheets, and vortex lines. The formulae obtained are the generalizations of the 


formulae of Helmholtz and Stokes for incompressible flow. The perturbation velocity 
s discontinuous at source layers and vortex sheets, and the discontinuities are 
xpressed in terms of the source and vorticity densities and the normals to the 


surfaces of discontinuity. 


1. The linearized equations of steady flow 
Ir U is the constant velocity of the undisturbed stream, a is the acoustic 
velocity, Vv is the perturbation particle velocity, § is the vorticity, Q is the 
source density of volume, and I is the unit tensor or idemfactor, then the 
inearized equations of flow for an inviscid fluid can be written in the form 
ef. (1, 2)) 
VAV C. V.wW Q, Ww (I— UU /a?).v. (34)7 
This system of equations is of the form treated in Part I (1) and hence 
solutions of the system can be written down immediately from the general 
esults contained therein. 
\sin Part Lit is convenient to introduce rectangular cartesian coordinates 
y,2.in this case such that the v-axis is parallel to, and in the same direction 
is, the velocity U. Let i, j, k be unit vectors parallel to the axes, and let 
VW U|/a) be the Mach number of the undisturbed stream; then the 
tensor I— UU /a? can be written in dyadic form as 
| —.W*)ii+-jj+kk, (35) 
which, by reference to Part I, section 3, shows that the system (34) is 
The equations aré iumbered so as to follow on from Part I. 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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elliptic for a subsonic main stream (.V < 1), and is hyperbolic for a super- 
sonic main stream (/ > 1), as was expected. The equations (34) are not 
valid for the physical problem when M 1, nor when is large, so that 
the conditions of Part I on the tensor (35) are always satisfied for the flow 
problem. The subsonic and supersonic cases have to be treated separately, 
as shown in Part I. 


2. Subsonic flow 


Part I, section 4.1, to be 


YW = Pii+jj+kk. (36) 
Then the system (34) is the special case of that considered in Part I for 
which f=v, g=w, P=, (det¥)!=8. (37) 


[fin addition Rg = |(x—2x,)?+f*(y—y,)?+f?(z—z,)?|!, 
then, from (20) of Part I, a solution of (34) at the point R, in the interior 
of a volume V bounded by a surface S is 


lf l\ 4, B f[mav(R)]A(R—R,) 5, 
, R pe Ww R V 1S es oe - 1 1S + 
v( 1) ron | n.w(R) (7) . io | Rs T 
Ss S , 
lf l ,, B® f (R—R,)AC(R) ,,, 
4 ( Vv 7 4. =... il IV, (3! 
i= | “(R) (,) 4 , | R3 a+ 7 


: 
where n is the unit outward normal to S. 

From the results of the next two sections it will be seen that the surface 
integrals in (38) can be interpreted as the respective contributions of a source 
distribution of surface density —n.w anda vorticity distribution, or vortex 
sheet, of surface density —n, von S. These surface distributions represent 
the effects of the sources and vortices outside S. However, as has been 
shown in Part I, section 4.3, nA v and n.w cannot be specified indepen- 
dently on S. If the solution is required for properly specified boundary 
values on S, it can be obtained by using suitable vector Green’s functions 
for S, as in Part I, section 4.4. 

The two volume integrals in (38) constitute a particular integral of the 
system (34), and they can be utilized for the determination of the velocity 
field due to source and vortex distributions. 

2.2. Inthe case when = 0 everywhere, and the only disturbance is that 
due to sources, the perturbation velocity field is given by 

i F ‘ . 
v(R,) = — Q(R)V( Re!) dV, (39) 
where the integration is over all space for which Q + 0. It is assumed that 
Q is such that the volume integral in (39) converges. 


2.1. In the subsonic case, put 8 = (1— M?)!, and take the tensor of 
I K 
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If the source distribution is confined to a layer, S, of vanishingly small 
thickness, in which Q is such that the surface density of sources, q, is finite, 
then the volume integral in (39) reduces to the surface integral 

v(R,) = — | q(R)V( Rg) dS, (40) 
a 
which gives the perturbation velocity field due to the source layer. 

[t is not difficult to show, either from (40), or directly from (34), that v 
and w are discontinuous on S. If v,, v, are the values of v on opposite sides 
of S at the same point, and n is the normal to S at this point, directed from 
side 1 to side 2, then 


nA (VvV.—V;) = 0, n.(W,—W),) = q, (41) 
—_" qn — q¥.n 
hence v,—V ; Ww,—W ; (42 
iad 5 . n.¥.n ‘ . a.F.2 


which gives the values of the discontinuities in terms of the source density, q, 


and the normal to the surface. Since 
n.¥. n 4 B? 


these discontinuities are finite, provided that q is bounded. 
In the physical flow problem, the layer is usually parallel to U at every 
point, in which case 


n.U Q), n.¥.n = 1, (43) 
and (42) gives the simple result, 


V.—V; w.—W, qn. (44) 


2.3. The perturbation velocity field due to a distribution of vorticity in 


the absence of sources is given by (38) as 


32 ¢ (R—R,) AG(R) 


v(R,) 
itr. Ri, 


dv, (45) 
where the integration is over all space for which $ 4 0. When M = 0, this 
reduces to Helmholtz’s classical result for the velocity field due to vorticity 
in an incompressible fluid, and represents its generalization for linearized 
subsonic compressible flow. 

When the vorticity is confined to a vanishingly thin layer, S, and is such 


that the surface density of vort icity is the finite vector w, (45) gives 


Q2 . (R —R,) A w(R) 79 
47r | Rs 


S 


v(R,) (46) 


Since the vortex lines must be continuous, w must be tangential to S at 
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every point of S. The velocity is discontinuous at such a vortex layer, or 
vortex sheet, as can be shown from (46), or from (34); if v,, V., and n have 
the meanings given to them in section 2.2, then 


Vv,—V, = nAw, (47) 


showing that only tangential components of v are discontinuous, and hence. 
since n.w = 0, Ps ; 

w@ = N/A (Vo—V)). (48) 
These results are independent of 8, and are the same as those for incom- 
pressible flow. 


Similarly, the perturbation velocity field for a vortex line, L, of strength 


wie BK f (R—R,) Ads 


wR) 4a R35 ; 


(49) 


te 


where ds is an element of L in the same direction as the vorticity. The 
integration in (49) is along L, and it can be shown that the line integral 
converges except on L. This formula is the generalization of that due to 
Stokes for an incompressible fluid. 


3. Supersonic flow 

3.1. In the case of supersonic flow, put B (M?—1)#, and take the 

tensor Y of Part I, section 5.1, to be 
WY — Bii—jj—kk. (50) 
Then the system (34) is the special case of that considered in Part I for which 
f = Vv, ¢ Ww, Pp «a (det PW)! B. (51) 
Also, let Ry = |(«—2,)?— B(y—y,)?— B?(z—z,)?}!. (52) 
In the physical flow problem, disturbances cannot be propagated upstream, 
and the effect of any disturbance located at the point R,, say, is confined 
to the domain Pe Pe 
, Ri; > 9, t> Xs; (53) 
which constitutes the domain of influence of R,. It follows that disturbances 


which affect the point R, must occur only at points that lie in the domain 


Ri, > 0, eB, (54) 


which constitutes the domain of dependence of R,; this domain is the 
interior of the characteristic half-cone, I, 


= : = 
R, Q, t< %,; (do) 


with vertex at R,. Thus, the boundary value problem considered in Part 5; 
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section 5.2, is the appropriate one for this case, and from (33) of Part I, 
(51), and (52), a solution for the perturbation velocity at the point R, is 





1 *° sep pee , _R | 
v(R,) = —5— | n-WiR)V(5 ) as [n wa (RR) gg) 
oe LR oT . : 
Ss re 
| * { | is B2* mi R ) P ¢(R) " ” 
= | ORV (z] 4V—— | 3, dV; (56) 


here S is that part of a surface on which n/ v and n.w are given which 
lies inside the characteristic half-cone [, V is the volume bounded by S 
and [’, and the star denotes that the finite part of the divergent integral 
is to be taken. 

Provided that the unit normal n to S satisfies the inequality 

a.¥.a > 6, (57) 

nv and n.w can be specified independently on S; otherwise this is not 
possible. 

As in the case of subsonic flow, the volume integrals in (56) give a par- 
ticular integral of the system (34), and they can be adapted to deal with 


cases of concentrated sources or vorticity. 


3.2. When the only disturbances in the whole field are due to sources, 


(56) gives 


} -F l : ich 
viRy) = 5. | QR)V( 2} dV, (58) 


the integration being over that part of the interior of [ for which Q + 0. 
For a source layer, S, of finite surface density q, the volume integral 
(58) reduces to a surface integral, and the perturbation velocity field is 
(R,) — -"f gRyv(2) as 59 
V(K,) = q( ) as, (59) 
- Ry 


. 
S 


the integration being over those parts of S inside’. As in the subsonic case, 
v and w are discontinuous at a source layer. The expressions for the 


discontinuities are formally identical with those given in (42), namely 


qn q¥.n 


: : (60) 
n.¥.n n.Y.n 


but here with the value for ¥ given by (50). In supersonic flow, n.¥.n = 0 
on characteristic surfaces, which shows that source layers cannot occur 
there, and that in the absence of source layers, the normal component of v 
and the tangential components of w can be discontinuous on, and only on, 
characteristic surfaces. 
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A result given by Robinson (3) for the discontinuity in v at a source layer 
agrees with (60) in the case when n. U = 0 (in which case (60) reduces to 
(44), as for subsonic flow), but appears to be incorrect when n.U + 0. 


3.3. When the only disturbance is due to a vorticity distribution, the 
velocity fields due to a volume distribution of vorticity, €, a vortex sheet, §, 
of surface vorticity w, and a vortex line, L, of strength K, can be derived 
as for the subsonic case, and are respectively 








2*;P - ; 
wy — BEF RY gy a 
— # B 
2* > P 
v(R,) = — .. = a oR) ag, (62) 
an 
B?K * — (R—R,) Ads 
v(R,;) = — — aps ? ma 
-— B 





where the integrations are over those parts of space, S, or L that lie inside 
the characteristic half-cone T. 

These generalizations of the results of Helmholtz and Stokes to the case 
of linearized supersonic flow have been obtained by Robinson (3), using a 
method based on a vector potential (the factors — B?/2z7 are omitted in 
error in ref. 3), and have been applied by him to calculate the downwash 
in the wake of a supersonic wing. 

The relations (47) and (48) between the velocity discontinuity and the 
surface vorticity are easily shown to hold for a supersonic vortex sheet as 
well as for a subsonic vortex sheet. 
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THE APPLICATION OF RELAXATION METHODS 
TO THE SOLUTION OF NON-ELLIPTIC PARTIAL 
DIFFERENTIAL EQUATIONS 


Il. THE SOLIDIFICATION OF LIQUIDS 


By D. N. DE G. ALLEN and R. T. SEVERN 
(Imperial College, London) 


Leceived 18 September 1951] 


SUMMARY 


In Part I of this series a method was described whereby the heat-conduction 


equation, . a 


could be attacked by relaxation in order to determine particular numerical solutions. 
The method is extended in this paper to determine the temperature distribution in 
a linear flow of heat in a liquid which solidifies as it cools through a given tempera- 
ture of solidification. The additional feature in this problem is the allowance that 
has to be made for the existence of latent heat of solidification, and particular 
attention must be paid to the determination of the progress of the solidification front 


is a function of time. 


1, ARELAXATIONAL method of solution of the heat-conduction equation, 


? Cn 
a aioe (1) 
ot Ox* 


was described in Part I of this series (1), where it was demonstrated that 
the temperature, v, in a linear flow of heat could be obtained numerically 
as a function of time, ¢, and distance, 2; in equation (1) « is the (uniform) 
diffusivityt of the medium. In this paper we extend the method so that 
it becomes capable also of giving the temperature in a linear flow of heat 
ina medium when a change of state occurs during that flow with a conse- 
quent liberation (or absorption) of latent heat; the method which will be 
described is illustrated by its application to determine the temperature in 

The diffusivity could be a function both of position (7) and of temperature (v); if a 
variation in either respect is taken into account, the form of the governing equation (1) is 
hanged and the computational problem becomes considerably more complicated. The 
mplications of a variation of the diffusivity, with respect to temperature in particular, have 

n extensively discussed by Eyres et al. (2) in a paper which describes the numerical 
solution of many heat-conduction problems by use of the differential analyser. Here we 
make no attempt to deal with other than a constant diffusivity, though allowance for its 


variation clearly merits further investigation. 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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a medium which extends from x = 0 tox = L, which is initially (at t = 0) 
liquid at its temperature of solidification V, where the end x = 0 is thermally 
insulated, and where at time ¢ = 0 the other end a = L is dropped to, and 
subsequently maintained at, a fixed temperature which for our purpose may 
be taken to be zero. As the medium cools it gradually solidifies from the 
end x = L and the progress of this solidification is seen in the (a, t)-plane 
as the locus of the point of solidification given by an equation of the form 

g = 2G). (2) 

The numerical solution to this problem is given in Fig. 1, which shows 
values of the temperature, v, at the nodes of a rectangular relaxation net 
in the (x, t)-plane; the net shown in Fig. 1 is the finest that was used in the 
computation of the solution, the mesh-length in the x-direction being 
h L/6 and in the t-direction H = L?/36x.2. The value taken for V was 45 
and that? for the latent heat of solidification (per unit volume) s 70-26 pe, 
where p is the density of the medium and ¢ is its specific heat; the total time- 
period covered by the solution is 77 = 5L?/3« 2. 

2. The relaxation method of solution of problems involving linear flows 
of heat is particularly suited to cases in which there is any generation of heat 
within the medium. In general we denote by A(z, t) the rate of generation of 
heat per unit volume at any point (and at any time) in the medium; the 
differential equation which then governs the temperature is 

ov Hy A , 
=« { ; (0) 
ot Cx* " pe 
and replaces (1), which holds only in those cases when there is no internal 
generation of heat. 

Now it was explained in Part I that equation (1) cannot be solved directly 
by relaxation methods for the reason that it is an equation of ‘marching’ 
type in the time coordinate; this difficulty was shown to be overcome by 
a preliminary transformation of (1) using the substitution 

» — w ae (4) 
ot Ox? 
to change the dependent variable from v to w. The same arguments hold in 
respect of (3) and it also must be transformed before a relaxational solution 
can be attempted. The same substitution (4) is made and so (3) is replaced 
7 @w_ate A _ 9, 5 
ot? éx* pe 


+ These particular numerical values were chosen in order that a comparison of our results 
could be made with those found from an exact analytical solution given by Jones (3). It 
will be seen in section 5 that it is not necessary to know s, p, and ¢ individually, but only 
to specify s/pc. 
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} When an appropriate solution of (5) is found for w, then v may be recovered 


results bv substitution back into (4) 
(3). I i ; 
aa 3. We use a rectangular relaxation net of sides h and H in the x- and t- 


rections respectively and, when a typical group of nodes of such a net are 
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denoted as shown in Fig. 2, the finite-difference approximation to (5) jg 
seen in the usual manner to be 


WetWy—2uy  K2(4w,+4w,—6wy—wy—w,,) Ag 
H? h4 pe 


which, after multiplication through by 2H?2, and with H chosent to be equal 
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0 








i 2 
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{ —e— A —i> 
oo —~e —> & 
11 3 0 1 d 
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Fic. 2 
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to h?/« v2, leads to the definition of a residual at a typical node © by the 
formula 

hy2 ; 
y—U'>,, —(A, hH) — . (6 
Kk 


ay Days Days — ’ ’ 
F, 2w,-+ 2w,+4w, + 4w,— 10w)—u 


where K (= xpc) is the conductivity of the medium. 

4. We next consider the physical interpretation of the factor 4,/H in 
the last term of the residual-formula (6). Fig. 3 represents a part of the 
relaxation net with 0 as a typical node; in cases when there is a generation 
of heat distributed in both time and distance, A, is the value of the rate ol 


+ The reason for the choice of this particular relation between h and H is give! 
Part I. 
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generation per unit volume at the node 0. But in cases when the generation 
of heat is concentrated either in time or position (spatially) then A,hH 
represents the total heat generated per unit cross-sectional area in the space- 
time interval shown as a shaded rectangle in Fig. 3. In the special case of the 
example which is treated numerically in this paper, the heat generated is 
itent heat and, if the locus (2) of the point of solidification passes through 
the shaded rectangle as shown also in Fig. 3, then the total heat generated 
within that rectangle’ is sd,, where d, is the projection on the x-axis of that 
part of the locus (2) which lies within the rectangle. Hence in such a case 
the formula (6) may be rewritten as 


j Shv2 
9 “11—% . ( 


K 


5. It is convenient to work numerically not in terms of w but in terms 


4a 4wu 


‘o y+ ZW, , LOwy—u 


~I 
~— 


3 


fa multiple cw/h? (= w’, say) (ef. Part I). For then the expression (7) is 
replaced by 


d,. sv2 
,Y Dan! Das’ / f m4 YY J 0 
I 2u,+ 2w,+4w, + 4u,— 10wp—wy— Uv}, 


h pe’ (8) 
nd in (8), in the last term, the factor d/h is a purely numerical quantity 
measuring the projection on the a-axis of the part of the solidification-locus 
shown in Fig. 3 as a fraction of the mesh-length in the 2-direction; similarly 
the factor (s\2)/(pc) is also a non-dimensionalt quantity whose value is a 
property of the medium only: in the numerical example worked in this paper 
the value of (sv2)/(pc) was taken to be 99-35 (ef. section 1). 

When values of w’ have been obtained at nodes of the relaxation net, 
corresponding values of the temperature v are immediately recovered from 
the finite-difference approximation to (4) 

Vg = wi + wy — 20s +s (wy—U). (9) 
VG 

6. In order that we should be able to solve equation (5) by relaxation to 
obtain a particular solution for w, certain boundary conditions must be 
stated; for the boundary of the rectangular region shown in Fig. | these 
onditions must comprise one on each of the sides = 0 and t T and two 

neach of the sides a 0 and a L. Some of these conditions have in 


effect been stated, in terms of v, in section 1; they are that 
on ft 0 v V ( "1 
on 2 0 CV Cx QO . 
, ink aca 
2)/(pe) is strictly non-dimensional only in the sense that it is of zero dimension in 


gth and time particular) and in mass; it is in fact a temperature. 
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or, in terms of w. 


ow Cw Bs 
on f¢ 0. — + 1 —— 45 | 
ct Ox" 
a3 | 
ou oPw | 
on 2 0), 1 4 : wy. (10 
cxct ox? | 
ou c*u 
and onz = L, — + «——. = 0 | 
ot Ox" 


The remaining three boundary conditions are arbitrarily chosen as simply 
as possible (cf. Part I) to be that 

on f¢ Tl, w 0 

onzxz=0. dw/ex oO 3. ‘VI 
and on x iY w= 6 
As a consequence of (11) being imposed. the original conditions (10) simplify 


to become that 


9 
ow Cw . 
on ft 0, 1 4 45 | 
‘ ‘ » 
ct Ox" | 
ou | 
on & 0, - 0 a (12) 
ox? 
Cw 
and on x L, -= 0 
Bx? 


In terms of w’ the finite-difference replacements for the boundary-conditions 
(11) and (12) are that 


onz=0, w, ws, and wy = wy, 
on x L, wW, 0 and wi)+ws 0 
(13) 
! iy ! f Dogs’ \alD 5/9 
on f 0, Wo— Wy + (W + W3— 29) V2 452 
and on f :,. = 0) 


these relations (13) are just sufficient for the elimination of all the fictitious 
values which occur in the use of (8) in order to calculate residuals at all nodes 
of the relaxation net. 

7. There only remains to describe how the solidification-locus (2) may be 
found, as an integral part of the solution, and how the factor d,/h, in the last 
term of (8), should be assessed at each node of the net. A process of system- 
atic trial and error is adopted: if any guessed position of the locus (2) is 
drawn on an accurate diagram of the net, the corresponding values of d, h 
holding at each node may be directly measured; using these values we can 
relax to determine a solution in terms of w’-values which satisfies the con- 
ditions (13) and makes residuals, defined by (8), vanish at all nodes. From 
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these w’-values, by substitution in (9), we can calculate nodal values of the 


temperature v: if the guessed position of the solidification-locus is correct 


we should find, as is manifestly required, that the value of v is 45 at ail nodes 


which lie in the region between the line ¢ 


t=5T/6 4 


oO 


15 at all nodes in the region between the guessed curve and the line ¢ 


RELAXATION 
| 


Vv 
\ 
\ 
ES \ 
\ 
~»\ 
\ 
\ 
\ 
\ 
\ 
\ 
\ 
yo a — 
r 
Kia. 4. 


0 and the locus, and is less than 


3A 


If the position of the curve is incorrectly guessed. then we should find 


values which do not conform with the above requirements: but in that 


event such v-values as are obtained give some indication of how a modified 


ocus should be drawn which can hardly fail to be an improvement on the 


original guess. Thus. if the v-values anywhere exceed 45 it is clear that the 
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position of the solidification-locus implied too rapid a generation of (latent) 
heat and therefore the modified curve should show a displacement away 
from the line t 0; conversely if the v-values tend to fall below 45 in the 
region between the curve and the line ¢ = 0, then the position of the curve 
implied too slow a generation of heat and the modified curve should be 
drawn nearer to the line ¢ = 0. The evidence of a very few trials enables a 
solidification-locus to be drawn which is very nearly correct, and subsequent 
minor. and local, modifications soon define it as accurately as may be 
desired. 

8. The solution? shown in Fig. 1 was obtained in the manner described 
in section 7, and, as a first example to be treated in this way, it is of interest 
that the results can be compared with those obtained from Jones’s exact 
analysis. For the particular numerical case dealt with here Jones gives the 
equation to the solidification-locus as 

(~a—L)? = 1-07x«t. (14) 
Fig. 4 reproduces (toa somewhat smaller scale) the curve which was obtained 
by the relaxation method and given in Fig. 1; it also shows Jones’s curve 
(14). It can be seen that the relaxational solution is thoroughly adequate 
for practical purposes; it could be improved, if necessary, by means of an 
advance to a finer relaxation net. 

Jones points out that the similar problem for heat-flow, which is radial 
(with cylindrical symmetry) instead of linear, has as yet proved intractable 
by exact analysis—his work is in fact an attempt to find an approximation 
to the total time of solidification of a cylinder of liquid, initially at its melting 
temperature and solidifying inwards from its surface, with cylindrical sym- 
metry. The extension of the relaxation method to solve the heat-conduction 
equation, in cylindrical coordinates, with allowance when required for the 
generation of latent heat. will be described in a subsequent paper. 


were obtained to four significant figures ; the work was stopped as soon as the modifications 


| gives values of the temperature v only. In the actual relaxation values of 


to the solidification-locus required for further refinement became inappreciable. 
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THE CALCULATION OF THE EFFICIENCY OF HEAT 
REGENERATORS 


By D. N. pe G. ALLEN (Imperial College, London) 
Received 18 September 1951] 


SUMMARY 


The ealeulation of the efficiency (or thermal ratio) of a heat regenerator requires 
£ { 
pair of simultaneous first-order partial differential 


equations in order to determine the temperatures of the fluid and of the matrix 


throughout one cycle of operation. A method is proposed in this paper which is 
partly iterative; its effectiveness is illustrated by three typical numerical results 
and, in particular, the question of convergence is discussed. 


1. Ir is the purpose of this paper to discuss only the numerical solution 
of certain partial differential equations which arise in the calculation of the 
thermal ratio of a heat regenerator; this entirely mathematical part of a 
ger investigation was brought to the author’s attention by Dr. S. 
Smoleniec, who has reported elsewhere (1) on the engineering problems 
involved and on the significance and application of the results. It is here 
sufficient, therefore. to define only notation and to state the governing 
equations; the physical quantities involved are denoted as follows:t 
L the length of the regenerator. 
v = distance along the regenerator measured from an origin at one end 
(QO ’ a. 
D = the duration of a single ‘blow’ of fluid (between reversals of flow) 
through the regenerator. 
time during a blow (0 2 <= Dp. 
the matrix surface area per unit length along the regenerator (per 


unit face area). 


h — the surface heat transfer coefficient between the matrix and the 
fluid. 

(; = the fluid mass rate of flow (per unit face area). 

m — the matrix mass per unit length along the regenerator (per unit face 
area) 


the specific heat of the fluid (per unit mass). 
c = the specific heat of the matrix (per unit mass). 
T’ = temperature of the matrix (a function of x and z). 


t = temperature of the fluid (a function of x and z). 


This notation closely follows that used in ref. 1. 


{Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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We use non-dimensional space and time coordinates given by the 


relations: ' ' 
X (sh Gc, )a | 


2 (] 
and Z (sh/me)z J 
the corresponding non-dimensional length and duration are 


A (sh Ge )L | 


p 
and 5 (sh/mc)D J” : 

2. One cycle of operation of the regenerator consists of the blow through 
it of air in one direction followed by a reverse blow in the other direction 
of gas: the air at entry is cold and may conveniently be taken to be at a 
temperature t = 0; the gas at entry is hot and taken to be at temperature 
t 1000. In its passage through the regenerator the air is heated through 
heat transfer from the (hotter) matrix; and conversely in its passage the gas 
is cooled through heat transfer to the (colder) matrix. The cycle is repre- 


sented diagrammatically in Fig. 1, in which it will be noticed that we 
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measure X always from the end of entry of the flow. The equations which 


covern the temperatures are then known to be 


oT 
t=27' 
al | . 
f ; : 
and : : 
in ax 


The boundary conditions to be satisfied by the solution are (Fig. 1) that 
t 0on AC 
t 1000 on DB’ 
and that the values of 7’ on A’B’ and on AB are equal 
3. We work primarily with the matrix temperature by eliminating ¢ from 
the equations (3) to derive a single equation governing 7’: 
eT eT or 
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this equation must be satisfied within each of the regions ABDCA and 
()B'A'C of Fig. 1, and in it X is measured to the right in the first of these 


regions, to the left in the second. 


We obtain the solution of (5) in the form of values of 7’ at the nodes of a 


rectangular net of sides a and 6 respectively in the X- and Z-directions ; 


Fig. 2 illustrates a typical mesh of such a net and with typical nodes denoted 


is shown, the finite-difference approximation to the governing equation 


5) can be easily verified to be 














= | |) 2+-a—5b\,.. 2—a—6\,. . 
i | | T; r } a (6) 
2+a+)/) 2+a+b 9+a+b 

The manner in which this equation is used to 6 5 
btain the numerical solutions will be explained ft 
in section 6. zZ b 
4. When values of the matrix temperature, 7’, - | 

have been found the fluid temperature, t, can be 8 

{/<=—a-— 


valuated directly by substitution in the first of 
3); but it is not, however, necessary to find ¢ 


explicitly in order that the thermal ratio be 








letermined. For, by definition, the thermal ratio ~Y 
sthe average change in temperature of the fluid FiG. 2. 
s it passes through the regenerator, measured as a fraction of the tem 
perature difference between the air and gas at entry (in this case, 1000); i.e. 
the thermal ratio 7 is given by 

9) 


t(A, Z) dZ, 


LOOOd 
. 


culated for the duration of a blow of air, and by 
re) 
| P P ‘ 
t(A.Z)dZ 
LOOOSd 
0 
for the duration of a blow of gas. 


Using the first of (3) we can transform these expressions into 


1 | . - 
| 7T(A.8)—T(A,0)+ | 110. 2) a2. (7) 
LOO00 | _ 
| . 
= T'(A,5)—T(A,0)+ | (A. 2) az}. (3) 
LOOOSd d 


irom which forms » may be easily evaluated once 7' has been determined. 
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5. Three cases have been worked numerically with results which are now 
summarized in a tabular form, as below. 


fir flou Gas flow 
A r?) n A r) 1] 
Ist net ‘708 “472 
Case | 2nd net 5°0 | 1:00 | *709 5°0 | Teo | *713 
3rd net *7O9 “710 
Ist net “840 ‘560 
Case II 2nd net 16°3 | 4°64 | 887 | 25-6 | 7:26 | -892 
3rd net “goo 895 
Ist net “905 72! 
Case III 2nd net 193 | 2°92 | +927 | 15°3 | 4°75 | °749 
3rd net ‘Q21 7604 


For each case a solution was obtained on three successive nets. Each time 
the first, or coarsest, net used had five mesh-lengths in the length of the 
regenerator and one mesh-length in the duration of a blow (between reversals 
of flow); the second net used had ten and two such mesh-lengths respectively; 
and the third, and finest, net used had similarly twenty mesh-lengths in the 
length of the regenerator and four in the duration of a blow. 

In case I, A and 6 have the same values for both the air and gas flows. 
This is the simplest cireumstance and implies equal total flows of air and gas 
per blow at equal mass rates of flow. In case IT, the values of A and 6 differ 
for the flows of air and gas but the ratio 5/A is the same: this corresponds to 
equal total flows per blow, but different mass rates of flow for air and gas. 
In case ITT, the values of A and 6 and of 5/A differ for the flows of air and gas: 
this means that unequal total flows of air and of gas per blow are being con- 
sidered, and also that the mass rates of flow are different. The more general 
circumstances of case III do not provide any extra difficulty in the numerical 
work required in the solution. Cases I, IT, and III are examples of continual 
reversals described by Smoleniec respectively as ‘symmetrical’, ‘asym- 
metrical’, and ‘unbalanced’. 

In each case the thermal ratio, », was determined using the three succes- 
sive nets. The convergence of the values of 7 from net to net is apparent 
from the results shown in the table, and they are sufficient for an estima- 
tion to be made of the exact value in each case. Thus, in case I, 7 is 
converging to 0-710: in case II it may be estimated that » is converging 
to 0-900. In case III there are two values for 7, one evaluated for each of 
the flows of air and gas respectively. For the air-flow the estimated exact 
value for 7 is 0-922, while for the gas-flow it is 0-780. 

These estimated values of 7 were computed on the assumption that the 
same ratio of errors on successive nets is perpetuated as the nets are made 
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increasingly finer. Thus if 7,, y., and ny are values obtained for y on three 


successive nets, the estimate of the exact value of 7 is given by 


1— 7) N2— 7 


2-7 %3— 7 


6. The distribution of 7'-values over the region A BB’A’'A (Fig. 1) is 
found by a computational process making repeated use of the relation (6). 
The process is started by an initial guess of the values at nodal points on A B; 
it once consistent values may be found at points on AC, for there ¢ = 0 and 
therefore, by the first of (3), 

Le; 
C Z 


so that 7 7,0". (9) 


With values for 7' at all nodal points of AB and AC, the relation (6) may 
be used to calculate successively 7'-values at points along each mesh-line 
parallel to AB until CD is reached. And then consistent values may be 





found at points on DB’, for there ¢ 1000 and therefore, by the first of 
} (3) again, ae 
oT 
eZ 


é i 1000, 


so that T 1000+ (T),— 1000)¢ Z (10) 


In (10), of course, an origin of Z is implied at D whereas in (9) the origin 
| of Zwas taken at A. Now with values of 7 at all nodal points of DC and 
| DB’ the relation (6) can be again used to calculate successively 7'-values at 
points along each mesh-line parallel to DC, until B’A’ is reached. 

\t this stage values of 7’ have been found over the whole region A BB’ A’ A: 
they will comprise the required solution if and only if their values along A’ B’ 
reproduce exactly the initially guessed values along AB. If this is not the 
case then those initial values must be altered and the process repeated until 
such a condition is achieved. Further details of this computational process 
will be mentioned in section 7, where its convergence is discussed: the results 
ire illustrated here in Fig. 3 which shows, for one (III) of the three cases 
worked, the final values of 7’ on the finest net used. There is no advantage to 
be gained by giving corresponding solutions for the other two cases (I and IT) 
rin giving any solutions on the coarser nets. 

7. There are two quite separate points of convergence to be mentioned; 
the first deals with the convergence, on any particular net, of the computa- 


tional process to yield an approximate numerical result peculiar to that net. 
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The second point deals with the convergence of such approximate results 
on successively finer nets towards a limit—the exact result. 

The first of these questions of convergence, in effect, is only concerned 
with the possibility of finding 7'-values on AB such that, when the process 
described in section 6 is carried through, these same values are reproduced 
on A’B’. A purely iterative procedure could be started with any guessed 
values on A B, the corresponding values on A’ B’ being calculated, these then 
put back as holding on A B and the cycle repeated until a steady state was 

















329 380 4x 46! 522, S76 62) 664 70S 744 780 B14 846 876 903 928 950 969 984 964 $09 
49 301 353, 406 457. 507, 5s5 60! 645 667. 727. 765 800 833 664 894 92! 946 968 986 9 
169 219 272 326 380 433 464 534 se 626 669 710 748 784 se 8s! 862 912 94) 968 99 
90) 136 68. 242 297 353 407 460 su ss59 606 650 691 730 "7 602 635 667 699. 932. 9 
8 S6 103 156 2 268 325 36! 4%6 487 538 586 63! 673 713 750 785 66 848 8 e 
37 23) 154 2i7. 279 340 398 454 $07. $57. 604 649 69! 730 767 sol 833 8 689. 92 9 
% 153. 226 295 359 419 475. $28 $77. 624 668 709 747 783 817 e468 876 902 925 944 9 
's 246 320 386, 445 soo. s 599. 644 687 727 764 799. 632 662 690. 915 938 95 97: 8. 
329 8 43) 48) 529 $76 6 664 705 44 780 B14 846 876 903 928 95 969 64 994 9 
Fic. 3 


assumed. Whether such an iterative process is convergent or not is largely 
immaterial, because in any event it would be quite unnecessarily laborious. 
Instead, a semi-relaxational approach easily enables the correct 7'-values 
on AB to be found after only two or three trials, and would indeed do so 
even if the iterative procedure actually diverged. For it is very soon appar- 
ent what effects are produced in the values on A’ B’ by changes to the values 
on A Band these can then be made accordingly to correct any discrepancies. 

The nearer the first guess of the values on A B is to the values which they 
eventually attain, the quicker, of course, is a solution found; if only for this 
reason it is desirable to work on successively finer nets, beginning with one 
which has very few nodal-points; for such a solution, found quickly on a 
coarse net, provides a relatively good initial set of values for the next finer 
net. All the three examples reported in this paper were worked in this way 
on three successive nets. Obtaining solutions on successive nets is also 
extremely important in order that the accuracy of their approximation to 
the exact solution be estimated; the change from net to net can be seen at 
a glance, so that it is immediately apparent when a sufficiently fine net has 
been used for any required order of accuracy. Thus, in section 5 it was 
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possible on the basis of values found from three successive nets to give 7 in 
each case to an accuracy which is believed to be within } per cent. 
This work was mostly carried out some three years ago and publication 


delayed while Dr. Smoleniec completed his entire investigation. 
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RELAXATION METHODS. A THREE-DIMENSIONAL 
MECHANICAL ANALOG) 


3y D. C. VAUGHAN 


(Messrs. L. G. Mouchel and Partners, Ltd., 
38 Victoria Street, London, S.W. 1) 
[Received 17 July 1951] 
SUMMARY 

An analogy is established between the approximate mathematical solution to 
any three-dimensional finite difference equation, obtained by relaxation methods, 
and the deflexions, in an arbitrarily chosen direction, of the nodes of a correspond 
ing, suitably loaded and tensioned, weightless, elastic three-dimensional net. 

The analogy is valid for any form of space lattice, but three regular forms are of 
particular use. 


1. RELAXATION methods, including a mechanical analogy, have been 
developed by Southwell and his co-workers to give solutions to a number 
of boundary-value problems where the governing equations are linear 
the boundary conditions on a closed boundary are given and linear, and 
the solution is unique and free from singularity. 

The required function is represented as the small transverse displace- 
ment of a tensioned, weightless, inextensible membrane, transversely 
loaded, and subject to constraining forces. The required continuous 
function, membrane, and constraints are then replaced, approximately, 
by a regular tensioned net, loaded and constrained at its nodes, and the 
transverse displacement (or required function) at each node is found by a 
systematic relaxation of the constraints perpendicular to the net until 
they all approximate to zero. 

The representation of the problem as the displacement of an inextensible 
net prevents the application of this particular analogy to a three-dimen- 
sional problem, because such a three-dimensional net cannot deflect. No 


such limitation applies, however, to the purely mathematical process of 


relaxation. The solution of Poisson’s equation with fixed boundary values 
has recently been achieved by Allen and Dennis (1) on a cubical net. It is 
to be expected that solutions will be deveioped in due course for any three- 
dimensional problem which can be formulated as a series of simultaneous 
equations between the values of the required function at the nodes of any 
suitable space lattice. 

Xepresentation of some three-dimensional problems by a mechanical 
analogy has so far been possible only where symmetry permits of mathe- 
matical reduction to a plane problem. 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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In what follows it is shown that the mathematical relaxation process 
for three-dimensional problems is equivalent to the small deflexion in an 
rbitrarily chosen direction, of the nodes of a suitably loaded and tensioned 
weightless elastic three-dimensional net. An identity is first established 
between the deflexion of the conventional plane and inextensible net, and 
that of a similar plane and elastic net, deflected in its plane. On this is 
based the argument for the extension to three-dimensional problems. 


Time is not included as an independent variable. 


2. For a regular, plane, inextensible net which is weightless and uni- 
formly tensioned, it has been shown by Christopherson and Southwell (2) 


that for small transverse displacements 


F ; p (w)—Nwy], (1) 
where, using the later notation of Southwell’s Relaxation Methods in 
Theoretical Physics (3): 
F, = the total force exerted on the node 0 by the N strings radiating 
from it. 
ss the string tension 
a the string length 


S (w) = the sum of the deflexions (w) of the N nodes distant a from 0. 


oy the deflexion of the node 0. 

\ uniform tension in all strings presupposes that the problem is governed 
vy Poisson’s equation throughout the whole domain. 

In the treatment of more complex problems, or of problems governed 

Poisson’s equation, but involving refraction at an interface, it has 

een assumed that the string tension varies in an arbitrary manner from 
string to string, or in the two portions of a string crossing an interface. In 
msequence, planar forces have been postulated acting on the nodes or 
strings in such a manner as to maintain equilibrium without displacement 
{the nodes in the plane of the net. Such forces do not affect the trans- 
erse displacement of the nodes, but as a consequence of the varying string 
tensions, equation (1) must now be replaced by 

, m - om 
F, |S (Tw) wy > (T)). (2) 


| « 
j 


1 ' \ aN 
If at any node all the radiating strings are not of equal length, equation 


F. 7 uw w, ty 7 (3) 
— a ~~ \a 


a 


Z become S 
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and this is the general equation governing the small transverse deflexio; 


of the nodes of any inextensible weightless tensioned plane net. 


3. Consider now a plane regular net formed of weightless tensione; 
elastic strings, and loaded at its nodes by parallel forces acting in an 
chosen direction in the plane of the net. 








6+60 


In the plane figure above, let the adjacent nodes A, B of such an unloaded 
net, in equilibrium, move in parallel directions through distances wy, w 
respectively, to the points A,, B, owing to the imposition of external forces 
fo. f, applied at the nodes in the direction of motion. Let the angle which 
these forces make with the normal to the string AB be 6. 

Further forces in the plane of the diagram, and perpendicular to fy, | 
must be postulated, but as these do not affect the displacements of A and B 
in the chosen direction A—A,, B—B,, they need not be considered here. 

Let the initial tension in the string, of initial length a, be 7’, and after 
deflexion let the tension increase to 7,. Suppose that the string obeys 


Hooke’s law. Then (tension)/(strain) E, a constant for the string under 
consideration. Hence 
- ee ; 
T, = T+ —(wy—w)sin 0. (4 
a 


In all that follows the displacements (w#,)—w) are assumed to be small in 
comparison with a. Then 


50 


and for equilibrium 


Fo T, sin(@+-86)— T sin 6. (6 
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From (4), (5), and (6), and by neglecting terms involving the second 
power Of (wW)—w)/a, 


Th 


Jo = _ T cos?6 + E sin?6). 


Suppose that .V strings, of lengths a, ay,...,a, meet at A; 0,a, T, and E 
will, in general, be different for each string, and, in the notation of equations 


1). (2), and (3), 


“Yb (T' cos?6é E sin?é)— Ws > ! (T cos?6- FE sin?6) ‘ (7) 

amt mv 

aN aN ‘ 

Equation (7) may be identified with (3) if 7’ for each string in equation 

3) is made equal to (7' cos*6+- F sin*@) for the same string in equation (7). 
Hence the transverse displacenient of the nodes of an inextensible tensioned 
net can be made identical wich the coplanar displacements, in an arbi- 
trarily chosen direction, of the nodes of a similar elastic tensioned net, by 
scribing suitable properties to the strings of the latter. 


4. Since the deflexion of the nodes of the elastic net is not confined to a 
direction perpendicular to the strings, there is no difficulty in extending the 
model analogy to three dimensions, and equation (7) is applicable without 
iteration to a three-dimensional net. 

\lso, since equation (7) does not depend on all strings in a net being of 
equal length, nor on a constant number of strings meeting at each node, 
it is applicable to a net arranged in any three-dimensional form, regular or 
irregular. In practice, to simplify computation the form of the net chosen 
should be such as will give identical relaxation patterns at all internal nodes. 
This limits the choice to those nets whose nodes form a cubic, cubic body- 
centred, or cubic face-centred space lattice. Each node then has 6, 8, and 
l2 strings respectively, of equal length, radiating from it. 
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HEAT CONDUCTION IN SOLIDS AS AN EIGENVALUE 
PROBLEM 
By G. D. WASSERMANN (King’s College, Newcastle upon Tyne) 
[Received 11 September 1951] 
SUMMARY 
Problems of linear heat conduction in solids are formulated as eigenvalue problems, 
This permits the rapid solution of standard problems. By the use of recent boundar 
perturbation techniques it also makes accessible to solution problems which cannot 
be solved by the usual methods. 
1. Introduction 
[n this paper we shall show that the majority of linear problems in the theory 
of heat conduction in solids can be treated by a unified method by reformu- 
lating these problems as eigenvalue problems.+ Our work is a development 
of a method due to Walters (1), who established a Laplace-transform 
relationship between the time-dependent Green’s function of the heat- 
conduction problem and the time-independent Green’s function of an 
associated elliptic differential equation. The approach presented in this 
paper will also apply, with slight modifications, to other differential equa- 
tions of theoretical physics. 

The present method has not only the advantage of unified presentation: 
it avoids the direct, and often cumbersome, application of Laplace trans- 
forms. Moreover, by formulating the problem as an eigenvalue problem we 
are in a position to obtain approximate solutions for problems which involve 
complicated boundaries and boundary conditions and which, apparently, 
cannot be solved by the usual methods (2). In these cases we can directly 
apply recently developed techniques (cf. 3, 4, 5). 

Some of the formulae which we shall derive are already known, but the 
majority are not to be found in the literature of the subject. It must be borne 
in mind that our treatment is of necessity formal, since convergence proofs 
are not possible for such a general formulation. These should, if possible, 
be given separately in each case. Whenever we have applied our methods 
to cases which have already been solved by other means, we have obtained 
identical results. 

2. Formal expressions for the time-dependent Green’s function 

We proceed to derive formal expressions for the Green’s function 
V(r,r’,t—t’) of the heat-conduction equation 

V2V = (1/x) eV /et (r Ar’,t st), (2.1) 
+ We assume in the following: (1) that the thermal properties are independent of tempera- 


ture; (2) that the boundary conditions are independent of, or linear in, the temperature; 
(3) that no change of phase or state occurs in the temperature range considered. 


(Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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where r is the position vector of a field point, r’ that of a unit pulse, and 
tand?t’ denote the times which refer to the field and unit pulse respectively. 
\ecording to Walters (1), V can be expressed in terms of the time-indepen- 
dent Green’s function G(r, r’,A) which obeys the same boundary conditions 


s ! and which satisfies the differential equation 


V2G 2G (rsr’); (2.2) 

4 becomes singular at r r’ in such a manner that over the domain D of 
ts definition : 

| {V2G(r,r’,A)—A*G(r,r’,A)} dr, & (2.3) 


D 


Once G is known, V is obtained as the inverse Laplace-transform of G, viz. 


Vir,r’,t—t’) = (1/2aix) | G(r,r’,¥ «)exp[A(t—t’)| dA (t > t’). 
: (2.4) 
We shall at first assume, without loss of generality, that the domain /) is 
finite. We shall further assume that V satisfies a suitable homogeneous 
boundary condition on the surface S of D. Later we shall show that infinite 
domains can be dealt with in an analogous manner. D is generally composed 
of several sub-domains DP. (7 i see n) which contain conductors of differ- 
ent thermal capacities x; (7 Be es n). It is to be understood that « in 
2.1) takes its appropriate value in D.. At the boundary of two adjacent 
sub-domains conditions of continuity of flux and temperature must be 
satisfied by V (and G@). 
Lets (r) anda, (n |, 2,...) be the normalized eigenfunctions and eigen 


alues respectively (in D) of the equation 


V7, (0) +Az wb, (r) = 0, (2.5) 


n 


where the ys,,(r) obey the same boundary and continuity conditions as 1’. 
For the sake of simplicity we shall assume that the eigenvalues are non- 
degenerate. It is readily shown, by a development similar to that used by 
Sommerfeld (6, pp. 183-4), that for a source of unit strength we have 
formally 
G(r,r’,A) ¥ &,,(r ps, (0 )[A2+Az2 |. (2.6) 


n 0 
Provided that the series on the right-hand side is boundedly convergent, 


we can insert it into (2.4) and obtain, formally, 


to 
~I 


Vireo" 4—5") > us, (r)b,,(r’ )exp| — Ka? (t t’)| ie => F 2. fe 
This expression forms the basis of our further development. We shall now 


systematically develop a theory of heat conduction based on (2.7). 
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3. The solution of initial-value problems 
We shall consider the general problem of Carslaw and Jaeger (2), chap. 13 
Case I. If v denotes the temperature function which satisfies (2.1) every. 
where at time ¢ > 0, and is such that 
v=f(r,0) at time? 0, 
and v d(r.t) onS at timet > 0, 
where f and ¢ are prescribed functions, then Minnigerode’s solution gives 
t 


v(r, t) Vir.r’.t)f(r’,0)dr,-+n | dt’ | ge’) oe, 2’, t') dS, 
e 
S 





Cc 
e 
D 0 


provided V satisfies the boundary condition 
| O ons. (3.2 
Suppose the eigenfunctions y,,(r) in (2.7) also satisfy (3.2): then fron 


(2.7) we obtain formally (e.g. if ¢ is independent of 1) 


vo(r,t) = > ¢,(r) ( ys, (r’) fie’, 0) dr, |exp(—d? t)4 
n=0 . 


D 





: » 
- ¥ [1—exp(— Kd? t)](1/A2 eb, (10) | d(r’) yp, (0") dS,.. (3.3 
n=0 R cv 
The first series represents the solution for the special case in which ¢ = ¢ 
which is the case usually given in the literature (cf. 6, p. 181). 

Case II. (2, p. 293.) If there is radiation into a medium at temperature ¢ 
the solution (3.3) applies, save that » must now satisfy the boundary 
condition 


opr’ 
: a : hv(r’) on 8S, (3.4 
and hence y,,(r) must satisfy 
oy, (0 
Cn ) hes, (v )} ons, (3.4 
Ov 


where / is a physical constant. 

As an example of the present method we treat a problem dealt with bj 
Carslaw and Jaeger (2, p. 306) by a different method. We consider : 
cylindrical region 0 < r < a, with zero surface temperature and unit instan- 
taneous source at f¢ 0. 


The solution of the appropriate eigenvalue problem 


d r 


rdr 


| d ( ele) | 2 b,,(r) 0. (3.5 
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which satisfies the equation 


ws QO atr a (3.6) 


ind which is finite on the axis of the cylinder, is 


Iola, r) 


dD (3.7) 
J,(A,, @) 
where J, and J, denote Bessel functions and A,, is a root of 
Ji(A,, @) 0. (3.8) 


Hence. from (3.3), we obtain at once 

JA, 7)" 

SFA, a)| . 
0 (3.9) 


which agrees with the solution of Carslaw and Jaeger. 


r'f(r’)Jg(A,, 7’) dr’, 


u(r t (2/a2) S exp( | 
n 1 


The theory for inhomogeneous media can be formulated in a similar 


manner. Here v satisfies div(«Vv) ov ct. (3.10) 


Hence the Green’s function is given by (2.7) provided that y,,(r) is a solu- 


tion of the eigenvalue problem 
div («Ves ) } d2 os, 0 (3.11) 


subject to appropriate boundary conditions. 

So far we have assumed a finite domain. In the case of an infinite or semi- 
infinite domain we obtain a continuous spectrum of eigenvalues, and hence 
the summation in (2.7) has to be replaced by an integral. As an example we 
consider the conduction of heat in a one-dimensional semi-infinite solid, 


oceupving the region x 0, with a bounding plane x 0, such that v 0 
itx = 0. The proper eigenvalue equation in this case is 
ds 
- + Ais 0, (3.12) 
dx 
subject to y = 0 at x = 0; and the normalized solutions which obey the 
boundary conditions are (p(A) eigenvalue density function) 
by ™p(A) |  sin(Az). (3.13) 


Hence the Green’s function is 


(2/7) ( sin(Ax)sin(Ax’ jexp| —«A?(t—t’)] dA 
7 4(«t)! |fexp| —(x’ —a)?/4«t|—exp[—(a’+-2)?/4xt]}, (3.14) 


which is precisely the expression used in section 20 of Carslaw and Jaeger’s 


book (2). 
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4. Continuous creation of heat 
We can also deal with the more general problem, which considers the 
creation of heat at a point r at the rate A(r,?/) per unit volume. In this case 
v satisfies the differential equation 
V2v (1/«) év/et+-(1/x)A(r, t). (4.] 
We obtain (cf. 1, p. 72), if the conditions are those which led to (3.3) from 
Minnigerode’s solution, 


t 


v(r,t) Vir.r’.t)f(r’,0) dr, +« | dt’ | dir’) as. 
J : : Cv 
t 
(1/x) [ dt’ [ V(rjr't—tA(r't) dry. (42 

0 D 

Consequently we get (3.1) (or (3.3) if d is independent of ¢) with an 

additional term 
t 


(1/«) = wb (vr) ( ys, (r’) ( A(r’,t’ )exp| r(t - t’)| dt’ dr, . (4.3 


n 0 


dD 0 


5. Stationary and moving sources 


Suppose that heat is generated at the point r’ by a stationary point source 
such that at time ¢’ the strength of the source is d(t’). Provided that (2.7) is 


boundedly convergent, we find that the temperature at a field point r at 
time ¢ is given by 


t t 


Vi(r,r’,t) [ V(r.r’,t—t’) dt’ > d,(r)y,(r’) [ f(t’ )exp| —KA2 (t—t’)| di 
0 n=0 0 (5.1 
In particular, for a source which is of constant strength q for ¢ > 0 we obtain 
V(r.r’,t) = (q/k) > &, (rh, (r’)A, >| 1—exp(— «A? t)]. (5.2 

n=0 


In many problems it is important to determine the effect of a ‘source 


aggregate’ such as line or band sources. If do, denotes a line or area element 


of sources, as the case may be, we obtain from (2.7) 
: t 
V (r,t) s u(r) ys, (r’) do, A(t’ )exp| KA? (t t’)| dt’. 5.3 
0 . . 


n 
o 0 


We can now make a transition to a moving point source. Suppose r’ moves 


so that r’(t’) is a known function of t’. In this case we obtain from (2.7) 0! 
the assumption of bounded convergence that 
: : 

Vir.t) = | V[r.r'(t’):t—t'] dt = > o,(r) | o,[r'(’yexp[— «3 (tt) dt. 

0 n 0 0 (5.4 
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The effect of an aggregate is similarly found to give 
e P 


Vir.t)= > y,(r) { [ ¢,[r'(e)Jexp[—d2(t—t’)] dt'do,. (5.5) 


From these formulae it is easy to deduce the special cases dealt with by 


Jaeger (7, 8 


The author wishes to thank Dr. J. D. Weston for reading this paper and 


for useful critical remarks. 
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SUMMARY 

The behaviour of Hamilton’s principle in non-holonomic systems has been investi- 
gated in two ways. H. Hertz (1), adopting generally accepted methods of the caleulus 
of variations, concluded that though there may be natural paths which satisfy th 
principle, there is an infinity of paths generated by the principle which are not natural, 
Hélder (2), on the other hand, using a method which is incompatible with the generally 
accepted formalism of the calculus of variations, was led to the conclusion that Hamil. 
ton’s principle is satisfied in non-holonomic systems. 

Hertz was a little unfortunate in his choice of an example; in the only instance ofa 
non-holonomic system to which he refers, all the natural paths do in fact satisfy 
Hamilton’s principle. This may have contributed to the frequent acceptance of 
Hélder’s incorrect theory, even today. In view of the misunderstandings which 
exist, it seems desirable to restate Hertz’s investigation, and to demonstrate that 
there are natural paths in non-holonomic systems which do not satisfy Hamilton’s 
principle; that is to say, natural paths for which f L dt has not a stationary value. 


1. Introduction 
HAMILTON’S principle, in non-holonomic systems, is expressed analytically 
by the requirement that (with the usual notation) 


1 
| L(q;,4,;,¢) dt (1) 


should be stationary with the conditions 


Qa. +8.) =o (Peer), @) 
k hey MS, 
(A repeated index, here and elsewhere, denotes summation.) 

The equations (2) are independent, and it is supposed that they do not 
admit of a complete set of integrals; for if they were completely integrable 
the system would be holonomic. 

The appropriate formalism for transforming this problem in the calculus 
of variations to a set of ordinary differential equations is well known. 
(For instance, Courant and Hilbert, Mathematische Physik, i, p. 191.) It 
appears to have been first established on a sound basis by A. Mayer (4). 
By it the equations (2) are to be multiplied by the Lagrange multipliers 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 4 (1952)] 
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\,() and added to the integrand; that is to say, the function L in condition 


1) is to be replaced by 
; 1,9 a ‘ 
L L(q;.9 »t) +A, (th QF (q;-t) T Q5(q;,8)} ( 1 .) ” 
" 2 


The customary process of variation, and partial integration, then gives 
73 £ 
1 


— Li) 8q; 0 (4) 


where the subscripts, g; and g;, denote partial derivatives. 

The quantities, 6g;, are independent: we have therefore, for a stationary 

value of the integral, 
ef nO Geet (5) 
dt 

These equations, with the r conditional equations (2), serve to determine 
the n coordinates, g;, and the r multipliers, A,,, as functions of the time. 

This method was adopted by H. Hertz (1) in his discussion of Hamilton’s 
principle in relation to non-holonomic systems, and it is the basis of the 
present inquiry 

Before continuing, an examination of the problem by Holder (2) and 
Voss (3) will be considered, so that the essential differences between the 
method which they employed and the normal formalism of the calculus 


of variations which has just been outlined may be revealed. 


2. Hélder’s method 

Whereas in the generally accepted method of the calculus of variations 
the comparison paths are required to satisfy the conditional equations, 
the displacements being free, in Hélder’s treatment this is reversed: the 
displacements are to satisfy the conditional equations, and it follows from 
the theory of Pfaff equations that the comparison paths do not satisfy 
them, a consequence which is indeed accepted by H6lder. The formalism 
is expressed by Voss in the following way. 

Varying | L dt without regard to the conditional equations, he derives 


| (z, se Sq, dt = 0. 


at 


. 


Combining this with the equations of condition Q* dq, 0 he deduces 
| i 4; 
the Lagrange equations 
_ dly, 
* dt 


L r,. a 


where the quantities A, are multipliers. 
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The last step is possible only if the displacements, 5q;, satisfy the con- 
ditional equations. 

No reasons, based on the theory of the calculus of variations, are given 
for the wide departures from the generally accepted formalism which are 
made. The method is briefly described here only because it appears stil] 
sometimes to be accepted as the correct approach to a consideration of 


Hamilton’s principle in relation to non-holonomic systems. 


3. The multiplicities of the family of paths derived from Hamilton’s 
principle, and the family of natural paths 
Accepting equations (5) as the differential equations generated by 
Hamilton’s principle, applied to a non-holonomic system, and substituting 
equations (3) in (5) we may express them in the form 


dL ea i 1,9 ie ) 
4% — DL +A, Q+A, Bi qg;+A, bh = 0 ( J (6) 
dt . ; : k eae r n 
where Bi; = Oa — av Bo = V:— Var (7) 


A subscript quantity, as distinct from an index, denotes partial differ- 
entiation for the quantity. 
The natural paths are known to be given by the solutions of Lagrange’s 


>quations ; 
ne dL;, k 2 2a n 

L,, py QF ; ' (8) 
dt k Deeks r<on 


where the quantities, j.,, are Lagrange multipliers. 

With each of the sets of equations, (6) and (8), we have to combine the 
conditional equations, (2). We have then, in each case, a set which suffices 
to determine the coordinates, g;, and the multipliers, A, or y;., as functions 
of the time. 

The multiplicities of the solutions of the two sets of equations, that is 
to say, the maximum number of independent constants of integration which 
can occur in them, will now be ascertained. Subsequently the effect on 
the relative multiplicities of the two solutions when the maximum number 
of constants is not realized, because there are ignorable coordinates present, 
will be considered. 

The set of equations (6) and (2) which embody the consequences of 
Hamilton’s principle contain derivatives, up to the second order, of the n 
coordinates, g;, and derivatives of the first order of the r multipliers, d,.. 
The maximum number of constants which can arise in the integrations is 
therefore 2n-+-r. They are the initial values of the coordinates and their 
velocities and the constants arising from the presence of the rates of change 
of the r multipliers. Since the initial coordinates and velocities must 
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satisfy the r conditional equations, (2), the maximum number of indepen- 
dent constants is 2x. This argument will hold, however, only if the con- 
ditional equations have no integral. For if there is an integral, then by 
algebraic manipulation we may arrange that one of the conditional equa- 


tions, say the mth, is a total differential equation; and if B}* and Bi‘ are 


the values of the quantities B*, and b‘,, defined by equations (7) when 


the conditional equations are rearranged in this way, the quantities 5” 
and Br" will be null for all values of iandj. The multiplier, A,,, will not then 


appear except in its rate of change, A,,: we are not therefore required to 
integrate this quantity, and one of the constants of integration will not 
arise. The full effect of the existence of integrals of the conditional equa- 
tions will be considered later: at the present stage we may state that the 
solutions of the equations derived from Hamilton’s principle, when the 
conditional equations have no integral, have at most 2x independent 
constants; that is to say, they are represented by a family of paths of 
multiplic ity o 

The 2” constants serve to define uniquely a path through any given pair 
of points. We observe also that there is an 00” multiplicity of paths through 
any point, in any direction compatible with the conditional equations. 
For if the » coordinates of the point are given, there is an n—r multiplicity 
of the initial velocities, since they are subject to the r equations of condition, 
2): the position and direction of the path element therefore absorb 2n—r 
of the 2n constants of integration, leaving r of them disposable. 

If there are integrals of the equations of condition, then for each integral 
the number of independent constants is reduced by 2. For if/(< r) integrals 
exist, say, we may, by algebraic manipulation, express the conditional 
equations as / equations in the form of total differentials, and r—l which 
have no integral. The/integralst of the first group may be used to eliminate 

of the coordinates: the equations (6) are thereby reduced in number to 

/, and r—/ non-integrable equations of condition remain. By the earlier 
rgument there are therefore 2(n—/) independent constants; that is to 
say, the number of constants is reduced by 2 for every integral of the con- 
ditional equations. 

Turning now to the Lagrange equations, the integration of equations 
8) and (2) introduces 2n constants, the number being less by r than in the 
integration of the equations derived from Hamilton’s principle because the 
tates of change of the r multipliers, 4;, do not occur. Again the number of 


independent constants is less by r, since the coordinates and their velocities 


\ constant must not be introduced on integration; for we suppose that the integral 


‘the expression of the geometrical relations in the system, and an arbitrary constant is 


therefore required 
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must satisfy the r equations of condition. There are therefore 2n—r inde- 
pendent constants and the multiplicity of the paths represented by the 
solutions is 00?"-". Through a given point, in a direction compatible with 
the conditional equations, there is clearly now one path only. 

If / integrals of the conditional equations exist, we may, by the earlier 
argument, reduce the problem to one in which there are n—/ coordinates 
and Lagrange equations, and r—/ non-integrable conditional equations, 
The maximum number of disposable constants is therefore 2(n—/)—(r—l), 
or 2n—r—l. The maximum number of constants of integration is con- 
sequently reduced by unity for each integral of the conditional equations, 

Collecting these results, we see that in a non-holonomic system having n 


coordinates whose rates of change are subject to 7 linear differential equa- 


tions, of which there are / (< r) integrals, the maximum multiplicity of 


the paths is 00%"~” by Hamilton’s principle, 0?”-"~ for the natural paths. 

The ratio of the maximum multiplicities of the Hamilton paths and the 
natural paths is therefore co"~’, 

If there are s ignorable coordinates in a system, since we can specify the 
paths of the system completely without introducing these coordinates 
explicitly, the number of constants is reduced by unity for each such 
coordinate. The multiplicity of the paths is now therefore «"~)-* by 
Hamilton’s principle, «0?”"~"~‘-* for the natural paths. The ratio of the 
multiplicities is unaffected. 


4. Natural motions which do not satisfy Hamilton’s principle 

We deduce from the relative multiplicities of the sets of paths of a system 
derived from Hamilton’s principle, and from Lagrange’s equations, that 
there are many paths in the former set which are not natural. We shall now 
inquire if there are natural paths which do not satisfy the principle: that 
is to say, which do not give { L dt a stationary value. 

Simple criteria by which those natural paths of a dynamical system which 
satisfy Hamilton’s principle may be distinguished from those that do not 
satisfy the principle, are easily derived when the system is a rigid body 
rolling on a rough horizontal plane under gravity. Here there are two 
non-integrable conditional equations giving the horizontal velocities, %, 
and #,, of the centre of mass, in terms of the Euler angles, 4, qs, 3, Say, 
specifying the orientation of the body in space. Observing that the time. 
and the coordinates of the centre of mass, do not enter explicitly, two of the 
equations (6), derived from Hamilton’s principle, are found to be 


mé,+A, = 0 (k = 1,2), (9) 


where m is the mass of the body. 
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The corresponding pair of Lagrange equations, (8), is 
mx, =p, (k 1, 2). (10) 
Substituting for #, and #, from the conditional equations, differentiated 
for the time, we see that X, and —,,, and X, and —pzy, are represented by the 
same functions of the Euler angles and their first and second derivatives. 
If, therefore, in the remaining three equations of the set (6), we substitute a 
solution of equations (8), they evidently reduce to 
i,j 1,2,3 
A, BE.g;+-A, Bk, = 0 J . (11) 
e= 1,3 
We may express the multipliers, A,, in terms of the Euler angles, their 
first derivatives, and two constants of integration, c, and c,, by integrating 
equations (9), and substituting for the velocities, #,., from the conditional 


equations. Omitting the terms A, B4, since the time does not enter explicitly, 


; : mT] ‘2.2 
AG; 9 » C;.) BF; q; 0 ( ie (12) 


we hav e 


as the required conditions that a solution of the Lagrange equations, 
q(t), that is to say a natural path, shall satisfy Hamilton’s principle. In 
attempting to satisfy the identities (12) the constants, c, and cy, are to be 
regarded as disposable. 

[fit is not possible to satisfy the identities, (12), by any choice of values 
of the constants, c, and c,, then the path does not satisfy Hamilton’s 
principle. For the expressions which are equated to zero in equations (6) are 
the factors of the variations, 4qg;, in the integrand of equation (4), and 
we have seen that along a natural path they reduce to A,(q;,q;.¢,) BK q;. 

1 
The integral | L dt therefore has a stationary value, for a natural path, 


in the class of non-holonomic systems considered, only if 


1 


, i,j 1,2,3 , 
[delay dys eu) BE, Gy git = 0 ("V5 ) (13) 


i 

Since the displacements, 5qg;, are arbitrary, the truth of the statement is 
evident. 

It will be observed that the integrand is the sum of the products of the 
multipliers, A,, and the bilinear covariant forms of the Pfaff equations 
which constitute the conditional equations. 

Taking, as a first example, a sphere rolling on a rough plane, for the 
natural paths the centre of mass moves in a straight line with constant 


velocity. Now integrating equations (9) we have 


MX, A, +-Cp, (k i, 2). 
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If we choose the constants of integration, c;,, to be equal to the constant 
linear momenta, m2;., we have A, = 0 and A, = 0. The conditions that the 
quantities, A, Bi,q,;, be identically zero along a natural path can therefore 
be satisfied by a suitable choice of the constants of integration, c, and c,, 
All natural paths therefore satisfy Hamilton’s principle. 


This was the only example referred to by Hertz in his discussion of 


Hamilton’s principle. The choice was not a happy one; for not only does 
the system provide no illustration of natural paths which do not satisf\ 
the principle, but in a sense it may be regarded as a holonomic system. 
For it is only at the instant of launch that forces, other than gravity, are 
acting at the point of contact of the sphere with the plane: throughout the 
subsequent path the natural motion is the same as that of a sphere launched 
in space under no forces, provided that the initial linear and angular 
velocities are chosen so that a point of the surface is instantaneously at 
rest. It may be that this unfortunate choice of an instance contributed to 
the failure of Hdlder, Voss, and others to accept Hertz’s conclusions. 

A system which provides an illustration of natural paths which do not 
satisfy Hamilton’s principle, as well as paths which do, is found in a thin 
circular disk rolling on a rough plane. The criteria (12) for this system are 
easily evaluated. Replacing the coordinates, q,, qs. 73 by 9, 4, %, where @ 
is the angle between the axis of the disk and the vertical, ¢ is the azimuth 
angle of the axis, and ¢ is the third Euler angle, the expressions on the 


left-hand sides of (12), are 


6-equation null 
-equation mary{aé sin @d—c, cos d—c, sin d} }. (14) 


y-equation mad{aésin 6—c, cosd—c,sin d} 


The expressions, (14), are identically zero if either 

(A) p= 0=¢4; 

(B) @ is constant, and c, and c, have zero values. 

There are known to be natural paths which satisfy the conditions (A), 
and the condition @ constant in (B). The former correspond to motions 
of the disk when it topples with the point of contact with the horizontal 
plane fixed, and the latter relate to rolling at constant inclination to the 
horizontal plane. These are therefore natural motions which satisfy 
Hamilton’s principle. 

On the other hand, there are natural motions of the disk which do not 
satisfy Hamilton’s principle, that is to say, motions represented by solutions 
of the Lagrange equations (8), which do not satisfy equations (6) derived 
from Hamilton’s principle, and which are not therefore consistent with a 


stationary value of | L dt. 
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Consider, for instance, the small oscillations of the disk about steady 
motion at constant inclination to the horizontal plane. In the steady motion 


let @= a, d = p, J+ C08 a n. In the oscillating motion let 


@ x4 € sin pt, 


where « is a small quantity. Then (Routh, Advanced Rigid Dynamics, 


9 


Section 2446. Ex. 3) gives the equations 
p?(k?-+ a?) p*tk?(1 +2 cos*ax)+-a? sin?! — np cos a(6k?+-a?)-4 
2n?(2k?+-a?)—gasina, (15) 


where k is the radius of gyration of the disk about a diameter. The values 


of u, x, and ” are required to satisfy the equation 
k?u? sin x cos «— (2k?+-a?)un sin a—ga cos « 0 (16) 
obtained by setting 6 = 0 in one of the equations of motion: they are 


therwise unrestricted 
We have seen that if the motion is to satisfy Hamilton’s principle the 
expressions (14) must be null. This requires, in this instance, that 
a@ sin @ c,cosd+c, sind. 
Substituting @ v-+esin pl, neglecting second-order terms, and choosing 


1 suitable origin for 4, this gives 


ep Sin « cos pt J(e}+c¢3)cos pt 
is the required condition, that is to say, 
p rv JF c3) eap sin & 
There are, in fact, natural motions satisfying these conditions. For if 
we set p = pw in equation (15), and eliminate the quantity n from (15) and 
16), we obtain an equation, quadratic in »?, which has real roots, of which 
ne is positive, for all values of the angle of inclination a. There is therefore 
a natural oscillatory motion, corresponding to each value of «, which 
satisfies Hamilton’s principle. If, however, we associate with any chosen 
value of «, a value of » which does not satisfy this quadratic equation, the 
condition p = p is not fulfilled, and the expressions (14) will be of the order 
«, as will therefore be the groups of terms, A, Bk,q;, which they represent. 
We see, in fact, that the equation (13) is, in the motion under consideration, 
1 
| mae(d d—y 54){ap sin a cos pt —,/(c}-+-c})cos pt} dt 0, 
j 
where (since only first-order terms in € are being retained) the angular 
velocities, 6 and ys, may be regarded as constant. For each value of the 
inclination of the disk there is therefore one natural path satisfying Hamil- 


ton’s principle, and oo! natural paths which do not. 
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5. Possible physical explanations of the failure of Hamilton’g 
principle 
Hertz (1) considered the possibility that the failure of Hamilton’g 
principle was due to some error in the mathematical formulation of rolling 
motion. He suggested that it might be basically incorrect to associat 


conservation of energy with rolling motion, which might necessarily 


involve some slipping. Other physical aspects could be examined in relation 
to the formulation: whether, for instance, the assumption of continuity; 
of rolling motion is justified. Little seems likely to be gained from such 
inquiries, however. The customary formulation, with whatever approximas 
tions are involved in it, does in fact, when associated with Lagrange’ 
method, give solutions which have a one-to-one correspondence with the 
natural paths. With Hamilton’s principle, on the other hand, it does not; 
The failure therefore arises in the principle itself, rather than in an im 
adequacy of the mathematical representation of the constraints. 
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